v0.13.1
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>
29
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(); }
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
51
64
71
74
77
78// using OpScaleL2 = MoFEM::OpScaleBaseBySpaceInverseOfMeasure<DomainEleOp>;
79
82#include <RigidBodies.hpp>
83#include <BasicFeTools.hpp>
84#include <MatrixFunction.hpp>
85
86// global variables
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>
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;
106
108
110
115
116private:
125
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
150
151 // SmartPetscObj<DM> dmElastic;
155
156public:
158 // global params
159 struct ProblemData {
160 int order;
161 PetscBool no_output;
162 PetscBool debug_info;
166 PetscBool scale_params;
167
170
171 PetscBool is_neohooke;
174 PetscBool is_ho_geometry;
175 PetscBool save_restarts;
176 PetscBool is_restart;
177
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
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
202
204 std::map<EntityHandle, int> mapBlockTets;
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";
241 };
242 };
243
245
246 struct MMonitor : public FEMethod {
247
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(
261 volSideFe->getOpPtrVector().push_back(
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(
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");
278 postProcFeSkin->getUserPolynomialBase() =
279 boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
280
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
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();
302 this, roller.methodOpForRollerPosition, roller.BodyDispScaled);
303 }
304 if (roller.methodOpForRollerDirection) {
305 roller.BodyDirectionScaled.clear();
307 this, roller.methodOpForRollerDirection,
308 roller.BodyDirectionScaled);
309 }
310 }
311 }
312
314 }
316
319
320 auto make_vtk = [&]() {
322 if (data.no_output || ts_step % data.output_freq != 0)
325 CHKERR postProcFe->writeFile("out_plastic_" +
326 boost::lexical_cast<std::string>(ts_step) +
327 ".h5m");
329 };
330
331 auto save_skeleton = [&]() {
333
335 CHKERR postProcFeSkeleton->writeFile(
336 "out_skeleton_" + boost::lexical_cast<std::string>(ts_step) +
337 ".h5m");
339 };
340
341 auto save_skin = [&]() {
343 if (data.no_output || ts_step % data.output_freq != 0)
345
347 CHKERR postProcFeSkin->writeFile(
348 "out_skin_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
350 };
351
352 double max_disp, min_disp;
353
354 auto print_max_min = [&](auto &tuple, const std::string msg) {
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);
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
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";
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
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);
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
465 CHKERR save_skeleton();
466
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
571 }
572
573 private:
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
590
593 CHKERR setUP();
595 CHKERR bC();
596 CHKERR OPs();
597 CHKERR tsSolve();
601}
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", "",
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", "",
628 PETSC_NULL);
629 CHKERR PetscOptionsBool(
630 "-is_fieldsplit_bc_only", "use fieldsplit for boundary only", "",
632 CHKERR PetscOptionsBool(
633 "-is_reduced_integration",
634 "use axiator split approach, with reduced integration", "",
636 CHKERR PetscOptionsBool("-test_user_base_tau", "test_user_base_for_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", "",
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);
683
684 char restart_file_name[255];
685 CHKERR PetscOptionsString("-restart_file", "restart file name", "",
686 "restart.dat", restart_file_name, 255,
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
698 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
699 "Neohooke material is not supported with plasticity.");
700
702}
703
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
791}
792
795 string block_name = "MAT_ELASTIC";
796 // TODO: implement blocks for plasticity
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
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)
828
829 cache = &mat_blocks.begin()->second;
830
832}
833
837
838 auto get_base = [this](bool has_hexes = false) {
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:
856 MOFEM_LOG("WORLD", Sev::inform)
857 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
858 break;
859 case DEMKOWICZ:
861 MOFEM_LOG("WORLD", Sev::inform)
862 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
863 break;
864 case BERNSTEIN:
866 MOFEM_LOG("WORLD", Sev::inform)
867 << "Set AINSWORTH_BERNSTEIN_BEZIER_BASE for displacements";
868 break;
869 default:
870 base = LASTBASE;
871 break;
872 }
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) {
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
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);
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) {
993 }
994 CHKERR simple->buildProblem();
996 CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
997
999}
1000
1003
1004 auto get_material_stiffness = [&]() {
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
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
1107 CHKERR commonDataPtr->calculateDiffDeviator();
1108
1109 CHKERR get_material_stiffness();
1110 if (data.is_neohooke)
1111 commonDataPtr->isNeoHook = true;
1112
1114}
1115
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;
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) {
1169 } else if (data.with_plasticity) {
1172 }
1173 // CHKERR DMoFEMPreProcessFiniteElements(dm, dirichlet_bc_ptr.get());
1174
1175 // forces and pressures on surface
1177 PETSC_NULL, "U");
1178 // nodal forces
1180 // edge forces
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(
1213
1215}
1216
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();
1259 bc_mng->getMergedBlocksMarker(vector<string>{"MOVE_", "ROTATE_"});
1260
1261 auto add_disp_bc_ops = [&]() {
1263 std::map<char, size_t> xyz_map = {{'X', 0}, {'Y', 1}, {'Z', 2}};
1264
1266 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
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
1338 };
1339
1340 auto add_ho_ops = [&](auto &pipeline) {
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());
1358 };
1359
1360 auto add_domain_base_ops = [&](auto &pipeline) {
1362 pipeline.push_back(
1364 if (data.is_ho_geometry)
1365 CHKERR add_ho_ops(pipeline);
1366
1367 pipeline.push_back(
1369 pipeline.push_back(
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) {
1384 "EP", commonDataPtr->plasticGradStrainPtr));
1385 pipeline.push_back(new OpCalculateScalarFieldGradient<3>(
1386 "TAU", commonDataPtr->plasticGradTauPtr));
1387 pipeline.push_back(
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
1398 };
1399
1400 auto add_domain_stress_ops = [&](auto &pipeline, auto mat_D_ptr) {
1402
1403 if (data.with_plasticity) {
1404 pipeline.push_back(new OpPlasticStress("U", commonDataPtr, mat_D_ptr));
1405
1406 pipeline.push_back(
1408
1409 } else
1410 pipeline.push_back(new OpStress("U", commonDataPtr, mat_D_ptr));
1411
1413 };
1414
1415 auto add_skeleton_base_ops = [&](auto &pipeline) {
1416 volSideFe = boost::make_shared<MoFEM::DomainSideEle>(mField);
1417
1418 volSideFe->getOpPtrVector().push_back(
1420 volSideFe->getOpPtrVector().push_back(
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) {
1429
1430 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
1432 pipeline.push_back(new OpSchurBegin("U", commonDataPtr));
1433
1434 if (data.is_large_strain) {
1435 if (data.is_neohooke)
1436 pipeline.push_back(
1438 else
1439 pipeline.push_back(
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) {
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
1464 };
1465
1466 auto add_domain_ops_rhs = [&](auto &pipeline, auto m_D_ptr) {
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(
1505 else
1506 pipeline.push_back(new OpCalculateLogStressTangent<false>(
1507 "U", commonDataPtr, m_D_ptr));
1508 pipeline.push_back(
1510 } else
1511 pipeline.push_back(
1513
1514 pipeline.push_back(new OpUnSetBc("U"));
1516 };
1517
1518 auto add_domain_ops_rhs_constraints = [&](auto &pipeline) {
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 }
1527 };
1528
1529 auto add_domain_ops_lhs_constraints = [&](auto &pipeline, auto m_D_ptr) {
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(
1541 pipeline.push_back(
1543 } else {
1544 pipeline.push_back(
1546 pipeline.push_back(
1548 }
1549
1550 pipeline.push_back(
1552 pipeline.push_back(
1554 pipeline.push_back(
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
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) {
1611 pipeline.push_back(
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(
1623 if (data.with_contact) {
1624 pipeline.push_back(
1626 pipeline.push_back(new OpCalculateHVecTensorTrace<3, BoundaryEleOp>(
1627 "SIGMA", commonDataPtr->contactTractionPtr));
1628 }
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(
1663
1664 pipeline.push_back(new OpRotatingFrameBoundaryLhs("U", "U", commonDataPtr,
1665 vol_side_fe_ptr));
1666 }
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));
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(
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
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
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) {
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
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(
1816 pipeline.push_back(new OpCalculateHVecTensorTrace<3, BoundaryEleOp>(
1817 "SIGMA", commonDataPtr->contactTractionPtr));
1818 pipeline.push_back(
1820 if (data.debug_info)
1821 pipeline.push_back(
1823 pipeline.push_back(new OpGetContactArea("SIGMA", commonDataPtr));
1824 data.update_roller->getRuleHook = integration_rule_bound;
1825 }
1826
1827 auto create_reactions_element = [&]() {
1829 Range &reaction_ents = commonDataPtr->reactionEnts;
1830
1831 auto get_reactions_meshset_ents = [&]() {
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 */) {
1854 "WORLD", Sev::warning,
1855 "Warning: Provided meshset (-reaction_id: %d) for calculating "
1856 "reactions is EMPTY! Calculating reactions for all FIX_ blocks.",
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
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(
1888 if (data.is_ho_geometry)
1889 CHKERR add_ho_ops(pipeline);
1890
1891 pipeline.push_back(
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(
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(
1917 else
1918 pipeline.push_back(new OpCalculateLogStressTangent<false>(
1919 "U", commonDataPtr, commonDataPtr->mtD));
1920 pipeline.push_back(
1922 } else
1923 pipeline.push_back(
1925
1927 };
1928
1930 CHKERR create_reactions_element();
1931
1933}
1934
1937
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) {
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);
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 = [&]() {
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 }
2004 };
2005
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);
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);
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,
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;
2093 };
2094
2095 auto postproc_dom = [&]() {
2097 CHKERR print_active_set(commonDataPtr->stateArrayPlast, "plastic",
2098 mField);
2100 };
2101 auto postproc_bound = [&]() {
2103 CHKERR print_active_set(commonDataPtr->stateArrayCont, "contact", mField);
2105 };
2106
2107 if (data.debug_info) {
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) {
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
2225 };
2226
2227 auto set_global_mat_and_vec = [&]() {
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);
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 = [&]() {
2249
2250 const MoFEM::Problem *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 }
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) {
2308 commonDataPtr->aoSEpEp = ep_data->getSmartRowMap();
2309 if (sigma_ao)
2310 commonDataPtr->aoSEpEp = sigma_ao;
2311
2312 commonDataPtr->blockMatContainerPtr =
2313 boost::make_shared<BlockMatContainer>();
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) {
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 = [&]() {
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);
2414 };
2415
2416 auto create_post_process_element = [&]() {
2418 postProcFe = boost::make_shared<PostProcVolumeOnRefinedMesh>(mField);
2419 postProcFe->generateReferenceElementMesh();
2421 postProcFe->getUserPolynomialBase() =
2422 boost::shared_ptr<BaseFunction>(new BubbleTauPolynomialBase());
2423 postProcFe->getOpPtrVector().push_back(
2425
2426 postProcFe->getOpPtrVector().push_back(
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(
2438 "EP", commonDataPtr->plasticStrainPtr));
2439 postProcFe->getOpPtrVector().push_back(
2441 postProcFe->getOpPtrVector().push_back(
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(
2451 else
2452 postProcFe->getOpPtrVector().push_back(
2454 commonDataPtr->mtD));
2455 }
2456 if (data.with_contact) {
2457 postProcFe->getOpPtrVector().push_back(
2459 "SIGMA", commonDataPtr->contactStressDivergencePtr));
2460 postProcFe->getOpPtrVector().push_back(
2462 "SIGMA", commonDataPtr->contactStressPtr));
2463
2464 postProcFe->getOpPtrVector().push_back(
2465 new OpPostProcContact("SIGMA", postProcFe->postProcMesh,
2466 postProcFe->mapGaussPts, commonDataPtr));
2467 }
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
2478 postProcFe->getOpPtrVector().push_back(
2479 new OpPostProcPlastic("U", postProcFe->postProcMesh,
2480 postProcFe->mapGaussPts, commonDataPtr));
2481
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");
2491 };
2492
2493 auto scatter_create = [&](auto coeff) {
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 = [&]() {
2509 boost::shared_ptr<MMonitor> monitor_ptr(
2512 boost::shared_ptr<ForcesAndSourcesCore> null;
2513 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
2514 monitor_ptr, null, null);
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
2577}
2578
2580
2582
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);
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("", "");
2638 } else
2639 CHKERR simple->loadFile("");
2640 }
2641
2642 // Add meshsets with material and boundary conditions
2643 CHKERR m_field.getInterface<MeshsetsManager>()->setMeshsetFromFile();
2644
2647 if (ex.data.scale_params)
2648 for (auto &mit : mat_blocks)
2649 CHKERR mit.second.scaleParameters();
2650
2651 CHKERR ex.runProblem();
2652 }
2654
2656}
static Index< 'p', 3 > p
std::string param_file
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.
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:304
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
Kronecker Delta class symmetric.
Kronecker Delta class.
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
FieldApproximationBase
approximation base
Definition: definitions.h:58
@ LASTBASE
Definition: definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ USER_BASE
user implemented approximation base
Definition: definitions.h:68
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
FieldSpace
approximation spaces
Definition: definitions.h:82
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ H1
continuous field
Definition: definitions.h:85
@ HDIV
field with continuous normal traction
Definition: definitions.h:87
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
auto integration_rule
constexpr auto t_kd
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: DMMMoFEM.cpp:747
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:373
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1144
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:533
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: DMMMoFEM.cpp:800
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
virtual const Field_multiIndex * get_fields() const =0
Get the fields object.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpSkeletonRhsPipeline()
Get the Op Skeleton Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual Field * get_field_structure(const std::string &name)=0
get field structure
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.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
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.
#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.
FTensor::Index< 'i', SPACE_DIM > i
double D
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
char mesh_file_name[255]
std::map< int, BlockParamData > mat_blocks
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixTensorTimesGradU< 3 > OpMiXLambdaGradURhs
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhs
static char help[]
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 3, 3 > OpMixUTimesLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixDivTimesVec< 3 > OpMixDivULhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixVecTimesDivLambda< 3 > OpMixUTimesDivLambdaRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpCentrifugalForce2
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpSpringRhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpMixDivTimesU< 3, 3, 3 > OpMixDivURhs
BlockParamData * cache
FormsIntegrators< BoundaryEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpSpringLhsNoFS
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, 3, 3, 1 > OpLogStrainMatrixLhs
FieldSpace space_test
constexpr bool TEST_H1_SPACE
EntityType zero_type
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, 3, 3, 0 > OpStiffnessMatrixLhs
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< USER_ASSEMBLE >::BiLinearForm< GAUSS >::OpMixTensorTimesGrad< 3 > OpLambdaGraULhs
Tensors class implemented by Walter Landry.
Definition: FTensor.hpp:51
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T * > &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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: DMMMoFEM.cpp:1003
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode set_contact_dm(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, SmartPetscObj< DM > &dm_contact)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
boost::weak_ptr< BlockMatContainer > block_mat_container_ptr
auto createTS(MPI_Comm comm)
auto createSmartGhostVector(MPI_Comm comm, PetscInt n, PetscInt N, PetscInt nghost, const PetscInt ghosts[])
Create smart ghost vector.
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
MoFEMErrorCode set_plastic_tau_dm(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, SmartPetscObj< DM > &dm_plastic_tau)
MoFEMErrorCode set_plastic_ep_dm(MoFEM::Interface &m_field, SmartPetscObj< DM > &dm, SmartPetscObj< DM > &dm_plastic)
Ddg< double, 3, 3 > tensor4_to_ddg(Tensor4< T, 3, 3, 3, 3 > &t)
struct OpSchurPreconditionMassContact< OpMassPrecHTensorBound > OpSchurPreconditionMassContactBound
struct OpSchurPreconditionMassContact< OpMassPrecHTensorVol > OpSchurPreconditionMassContactVol
constexpr AssemblyType A
Definition: plastic.cpp:35
constexpr IntegrationType G
Definition: plastic.cpp:36
constexpr auto field_name
static constexpr int approx_order
MoFEMErrorCode getOptionsFromCommandLine()
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode preProcess()
function is run at the beginning of loop
boost::shared_ptr< DomainSideEle > volSideFe
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
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)
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< OpContactTools::CommonData > commonDataPtr
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProcFe
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcFeSkin
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProcFeSkeleton
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
boost::shared_ptr< DomainEle > calc_reactions
boost::shared_ptr< DomainEle > contact_pipeline
std::map< EntityHandle, int > mapBlockTets
boost::shared_ptr< BoundaryEle > update_roller
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
FieldApproximationBase base
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
MoFEMErrorCode postProcess()
boost::ptr_vector< DataFromMove > bcData
MoFEMErrorCode loadMeshBlockSets()
boost::ptr_map< std::string, EdgeForce > edge_forces
boost::shared_ptr< OpContactTools::CommonData > commonDataPtr
MoFEM::Interface & mField
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
MoFEMErrorCode createCommonData()
boost::ptr_map< std::string, NeumannForcesSurface > neumann_forces
static MoFEMErrorCode getAnalysisParameters()
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProcFe
boost::shared_ptr< DomainSideEle > volSideFe
SmartPetscObj< DM > dmContact
static MoFEMErrorCode myMeshPartition(MoFEM::Interface &m_field)
boost::shared_ptr< DirichletDisplacementRemoveDofsBc > dirichlet_bc_ptr
static ProblemData data
ContactPlasticity(MoFEM::Interface &m_field)
boost::ptr_map< std::string, NodalForce > nodal_forces
MoFEMErrorCode checkResults()
MoFEMErrorCode runProblem()
SmartPetscObj< DM > dmEP
boost::shared_ptr< DomainEle > feLambdaLhs
SmartPetscObj< DM > dmTAU
boost::shared_ptr< DomainEle > feLambdaRhs
chrono::high_resolution_clock::time_point stop
chrono::high_resolution_clock::time_point start
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
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
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.
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.
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
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
static MoFEMErrorCode applyScale(const FEMethod *fe, boost::ptr_vector< MethodForForceScaling > &methods_op, VectorDouble &nf)
Simple interface for fast problem set-up.
Definition: BcManager.hpp:23
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:218
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:165
Managing BitRefLevels.
Managing BitRefLevels.
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
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
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:82
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
structure for User Loop Methods on finite elements
Provide data structure for (tensor) field approximation.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
Interface for managing meshsets containing materials and boundary conditions.
Calculate HO coordinates at gauss points.
Calculate divergence of tonsorial field using vectorial base.
Calculate tenor field using vectorial base, i.e. Hdiv/Hcurl.
Calculate trace of vector (Hdiv/Hcurl) space.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Get field gradients at integration pts for symmetric tensorial field rank 2.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Calculate normals at Gauss points of triangle element.
transform Hdiv base fluxes from reference element to physical triangle
Set indices on entities on finite element.
Apply contravariant (Piola) transfer to Hdiv space for HO geometry.
Set inverse jacobian to base functions.
transform local reference derivatives of shape function to global derivatives if higher order geometr...
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
keeps basic data about problem
boost::shared_ptr< SubProblemData > subProblemData
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:334
intrusive_ptr for managing petsc objects
TS ts
time solver
PetscReal ts_t
time
PetscInt ts_step
time step number
Vec & ts_u
state vector
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Common data]
Force scale operator for reading two columns.