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