v0.14.0
nonlinear_dynamics.cpp
Go to the documentation of this file.
1 /**
2  * \file nonlinear_dynamics.cpp
3  * \example nonlinear_dynamics.cpp
4  * \ingroup nonlinear_elastic_ele
5  *
6  * \brief Non-linear elastic dynamics.
7
8  \note For block solver is only for settings, some features are not implemented
9  for this part.
10
11  \note This is implementation where first order ODE is solved, and displacements
12  and velocities are independently approximated. User can set lowe approximation
13  order to velocities. However this method is inefficient comparing to method
14  using Alpha method sor second order ODEs. Look into tutorials to see how to
15  implement dynamic problem for TS type alpha2
16
17
18  */
19
20
21
22 #include <BasicFiniteElements.hpp>
23 using namespace MoFEM;
24
25 #include <ElasticMaterials.hpp>
27
28 static char help[] = "...\n\n";
29
30 struct MonitorPostProc : public FEMethod {
31
34  std::map<int, NonlinearElasticElement::BlockData> &setOfBlocks;
36  &feElasticEnergy; ///< calculate elastic energy
38  &feKineticEnergy; ///< calculate elastic energy
39
40  bool iNit;
41
42  int pRT;
43  int *step;
44
46  MoFEM::Interface &m_field,
47  std::map<int, NonlinearElasticElement::BlockData> &set_of_blocks,
48  NonlinearElasticElement::MyVolumeFE &fe_elastic_energy,
49  ConvectiveMassElement::MyVolumeFE &fe_kinetic_energy)
50  : FEMethod(), mField(m_field), postProc(m_field),
51  setOfBlocks(set_of_blocks), feElasticEnergy(fe_elastic_energy),
52  feKineticEnergy(fe_kinetic_energy), iNit(false) {
53
54  double def_t_val = 0;
55  const EntityHandle root_meshset = mField.get_moab().get_root_set();
56
57  Tag th_step;
58  rval = m_field.get_moab().tag_get_handle(
59  "_TsStep_", 1, MB_TYPE_INTEGER, th_step,
60  MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_t_val);
61  if (rval == MB_ALREADY_ALLOCATED) {
62  MOAB_THROW(m_field.get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
63  (const void **)&step));
64  } else {
65  MOAB_THROW(m_field.get_moab().tag_set_data(th_step, &root_meshset, 1,
66  &def_t_val));
67  MOAB_THROW(m_field.get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
68  (const void **)&step));
69  }
70
71  PetscBool flg = PETSC_TRUE;
72  CHK_THROW_MESSAGE(PetscOptionsGetInt(PETSC_NULL, PETSC_NULL,
73  "-my_output_prt", &pRT, &flg),
74  "Can not get option");
75  if (flg != PETSC_TRUE) {
76  pRT = 10;
77  }
78  }
79
82
83  if (!iNit) {
85  if(mField.check_field("MESH_NODE_POSITIONS"))
86  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", postProc, true, false, false,
87  false);
92
93  std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
94  setOfBlocks.begin();
95  for (; sit != setOfBlocks.end(); sit++) {
96  postProc.getOpPtrVector().push_back(new PostProcStress(
97  postProc.postProcMesh, postProc.mapGaussPts, "DISPLACEMENT",
98  sit->second, postProc.commonData, true));
99  }
100
101  iNit = true;
102  }
103
104  if ((*step) % pRT == 0) {
105  CHKERR mField.loop_finite_elements("DYNAMICS", "MASS_ELEMENT", postProc);
106  std::ostringstream sss;
107  sss << "out_values_" << (*step) << ".h5m";
108  CHKERR postProc.writeFile(sss.str().c_str());
109  }
110
111  feElasticEnergy.ts_ctx = TSMethod::CTX_TSNONE;
112  feElasticEnergy.snes_ctx = SnesMethod::CTX_SNESNONE;
113  CHKERR mField.loop_finite_elements("DYNAMICS", "ELASTIC", feElasticEnergy);
114
115  feKineticEnergy.ts_ctx = TSMethod::CTX_TSNONE;
116  CHKERR mField.loop_finite_elements("DYNAMICS", "MASS_ELEMENT",
117  feKineticEnergy);
118  double E = feElasticEnergy.eNergy;
119  double T = feKineticEnergy.eNergy;
120  MOFEM_LOG_C(
121  "DYNAMIC", Sev::inform,
122  "%d Time %3.2e Elastic energy %3.2e Kinetic Energy %3.2e Total %3.2e\n",
123  ts_step, ts_t, E, T, E + T);
124
126  }
127
131  }
132
136  }
137 };
138
139 struct MonitorRestart : public FEMethod {
140
141  double *time;
142  int *step;
144  int pRT;
145
146  MonitorRestart(MoFEM::Interface &m_field, TS ts) : mField(m_field) {
147  double def_t_val = 0;
148
149  const EntityHandle root_meshset = mField.get_moab().get_root_set();
150
151  Tag th_time;
152  rval = m_field.get_moab().tag_get_handle(
153  "_TsTime_", 1, MB_TYPE_DOUBLE, th_time,
154  MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_t_val);
155  if (rval == MB_ALREADY_ALLOCATED) {
156  rval = m_field.get_moab().tag_get_by_ptr(th_time, &root_meshset, 1,
157  (const void **)&time);
158  MOAB_THROW(rval);
159  ierr = TSSetTime(ts, *time);
160  CHKERRABORT(PETSC_COMM_WORLD, ierr);
161  } else {
162  rval = m_field.get_moab().tag_set_data(th_time, &root_meshset, 1,
163  &def_t_val);
164  MOAB_THROW(rval);
165  rval = m_field.get_moab().tag_get_by_ptr(th_time, &root_meshset, 1,
166  (const void **)&time);
167  MOAB_THROW(rval);
168  }
169  Tag th_step;
170  rval = m_field.get_moab().tag_get_handle(
171  "_TsStep_", 1, MB_TYPE_INTEGER, th_step,
172  MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_t_val);
173  if (rval == MB_ALREADY_ALLOCATED) {
174  rval = m_field.get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
175  (const void **)&step);
176  MOAB_THROW(rval);
177  } else {
178  rval = m_field.get_moab().tag_set_data(th_step, &root_meshset, 1,
179  &def_t_val);
180  MOAB_THROW(rval);
181  rval = m_field.get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
182  (const void **)&step);
183  MOAB_THROW(rval);
184  }
185
186  PetscBool flg = PETSC_TRUE;
187  ierr = PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_output_prt", &pRT,
188  &flg);
189  CHKERRABORT(PETSC_COMM_WORLD, ierr);
190  if (flg != PETSC_TRUE) {
191  pRT = 10;
192  }
193  }
194
197
198  (*time) = ts_t;
199  // if(pRT>0) {
200  // if((*step)%pRT==0) {
201  // std::ostringstream ss;
202  // ss << "restart_" << (*step) << ".h5m";
203  // CHKERR
204  // mField.get_moab().write_file(ss.str().c_str()/*,"MOAB","PARALLEL=WRITE_PART"*/);
205  //
206  // }
207  // }
208  (*step)++;
210  }
211
215  }
216
220  }
221 };
222
223 // See file users_modules/elasticity/TimeForceScale.hpp
224 #include <TimeForceScale.hpp>
225
226 int main(int argc, char *argv[]) {
227
228  const string default_options = "-ksp_type fgmres \n"
229  "-pc_type lu \n"
230  "-pc_factor_mat_solver_type mumps \n"
231  "-mat_mumps_icntl_20 0 \n"
232  "-ksp_atol 1e-10 \n"
233  "-ksp_rtol 1e-10 \n"
234  "-snes_monitor \n"
235  "-snes_type newtonls \n"
236  "-snes_linesearch_type basic \n"
237  "-snes_max_it 100 \n"
238  "-snes_atol 1e-7 \n"
239  "-snes_rtol 1e-7 \n"
240  "-ts_monitor \n"
241  "-ts_type alpha \n";
242
243  string param_file = "param_file.petsc";
244  if (!static_cast<bool>(ifstream(param_file))) {
245  std::ofstream file(param_file.c_str(), std::ios::ate);
246  if (file.is_open()) {
247  file << default_options;
248  file.close();
249  }
250  }
251
252  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
253
254  // Add logging channel for example
255  auto core_log = logging::core::get();
258  LogManager::setLog("DYNAMIC");
259  MOFEM_LOG_TAG("DYNAMIC", "dynamic");
260
261  try {
262
263  moab::Core mb_instance;
264  moab::Interface &moab = mb_instance;
265
266  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
267  auto moab_comm_wrap =
268  boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
269  if (pcomm == NULL)
270  pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
271
273  char mesh_file_name[255];
274  PetscBool is_partitioned = PETSC_FALSE;
275  PetscBool linear = PETSC_TRUE;
276  PetscInt disp_order = 1;
277  PetscInt vel_order = 1;
278  PetscBool is_solve_at_time_zero = PETSC_FALSE;
279
280  auto read_command_line_parameters = [&]() {
282  PetscBool flg = PETSC_TRUE;
283  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
284  mesh_file_name, 255, &flg);
285  if (flg != PETSC_TRUE)
286  SETERRQ(PETSC_COMM_SELF, 1, "Error -my_file (mesh file needed)");
287
288  // use this if your mesh is partitioned and you run code on parts,
289  // you can solve very big problems
290  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_is_partitioned",
291  &is_partitioned, &flg);
292
293  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-is_linear", &linear,
294  PETSC_NULL);
295
296  enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, LASBASETOP };
297  const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier"};
298  PetscInt choice_base_value = BERNSTEIN_BEZIER;
299  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
300  LASBASETOP, &choice_base_value, PETSC_NULL);
301  if (choice_base_value == LEGENDRE)
303  else if (choice_base_value == LOBATTO)
304  base = AINSWORTH_LOBATTO_BASE;
305  else if (choice_base_value == BERNSTEIN_BEZIER)
307
308  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_disp_order",
309  &disp_order, &flg);
310  if (flg != PETSC_TRUE)
311  disp_order = 1;
312
313  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_vel_order",
314  &vel_order, &flg);
315  if (flg != PETSC_TRUE)
316  vel_order = disp_order;
317
318  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL,
319  "-my_solve_at_time_zero",
320  &is_solve_at_time_zero, &flg);
321
323  };
324
325  auto read_mesh = [&]() {
327  if (is_partitioned == PETSC_TRUE) {
328  // Read mesh to MOAB
329  const char *option;
330  option = "PARALLEL=BCAST_DELETE;"
331  "PARALLEL_RESOLVE_SHARED_ENTS;"
332  "PARTITION=PARALLEL_PARTITION;";
334  } else {
335  const char *option;
336  option = "";
338  }
340  };
341
344
345  MoFEM::Core core(moab);
346  MoFEM::Interface &m_field = core;
347
348  // ref meshset ref level 0
349  BitRefLevel bit_level0;
350  bit_level0.set(0);
351  EntityHandle meshset_level0;
352  CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
353  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
354  0, 3, bit_level0);
355  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
356  bit_level0, BitRefLevel().set(), meshset_level0);
357
358  // Fields
360  3, MB_TAG_SPARSE, MF_ZERO);
362  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
363  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
364  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
365  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
366
367  bool check_if_spatial_field_exist = m_field.check_field("DISPLACEMENT");
368  CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_SPARSE,
369  MF_ZERO);
370  // add entities (by tets) to the field
372
373  // set app. order
374  CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", disp_order);
375  CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", disp_order);
376  CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", disp_order);
378  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", disp_order);
379  else
380  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
381
382  // Add nodal force element
386  // Add fluid pressure finite elements
387  FluidPressure fluid_pressure_fe(m_field);
390  false, false);
392  "DISPLACEMENT", PETSC_NULL, false, true);
393
394  // Velocity
395  CHKERR m_field.add_field("VELOCITY", H1, base, 3, MB_TAG_SPARSE, MF_ZERO);
397
398  CHKERR m_field.set_field_order(0, MBTET, "VELOCITY", vel_order);
399  CHKERR m_field.set_field_order(0, MBTRI, "VELOCITY", vel_order);
400  CHKERR m_field.set_field_order(0, MBEDGE, "VELOCITY", vel_order);
402  CHKERR m_field.set_field_order(0, MBVERTEX, "VELOCITY", vel_order);
403  else
404  CHKERR m_field.set_field_order(0, MBVERTEX, "VELOCITY", 1);
405
406  CHKERR m_field.add_field("DOT_DISPLACEMENT", H1, base, 3, MB_TAG_SPARSE,
407  MF_ZERO);
409  CHKERR m_field.set_field_order(0, MBTET, "DOT_DISPLACEMENT", disp_order);
410  CHKERR m_field.set_field_order(0, MBTRI, "DOT_DISPLACEMENT", disp_order);
411  CHKERR m_field.set_field_order(0, MBEDGE, "DOT_DISPLACEMENT", disp_order);
413  CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_DISPLACEMENT",
414  disp_order);
415  else
416  CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_DISPLACEMENT", 1);
417
418  CHKERR m_field.add_field("DOT_VELOCITY", H1, base, 3, MB_TAG_SPARSE,
419  MF_ZERO);
421  CHKERR m_field.set_field_order(0, MBTET, "DOT_VELOCITY", vel_order);
422  CHKERR m_field.set_field_order(0, MBTRI, "DOT_VELOCITY", vel_order);
423  CHKERR m_field.set_field_order(0, MBEDGE, "DOT_VELOCITY", vel_order);
425  CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_VELOCITY", disp_order);
426  else
427  CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_VELOCITY", 1);
428
429  // Set material model and mass element
430  NonlinearElasticElement elastic(m_field, 2);
431  ElasticMaterials elastic_materials(m_field);
432  CHKERR elastic_materials.setBlocks(elastic.setOfBlocks);
435  // NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<double>
436  // st_venant_kirchhoff_material_double; CHKERR
439  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeRhs(), true, false,
440  false, false);
441  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeLhs(), true, false,
442  false, false);
444  false, false, false);
445  CHKERR elastic.setOperators("DISPLACEMENT", "MESH_NODE_POSITIONS", false,
446  true);
447
448  // set mass element
449  ConvectiveMassElement inertia(m_field, 1);
450  // CHKERR inertia.setBlocks();
451  CHKERR elastic_materials.setBlocks(inertia.setOfBlocks);
453  "DISPLACEMENT");
455  "DISPLACEMENT");
456
458  {
459  string name = "-my_accelerogram";
460  char time_file_name[255];
461  PetscBool flg;
462  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, name.c_str(),
463  time_file_name, 255, &flg);
464  if (flg == PETSC_TRUE) {
465  inertia.methodsOp.push_back(new TimeAccelerogram(name));
466  }
467  }
468
469  // damper element
470  KelvinVoigtDamper damper(m_field);
471  CHKERR elastic_materials.setBlocks(damper.blockMaterialDataMap);
472  {
473  KelvinVoigtDamper::CommonData &common_data = damper.commonData;
474  common_data.spatialPositionName = "DISPLACEMENT";
475  common_data.spatialPositionNameDot = "DOT_DISPLACEMENT";
478  "DISPLACEMENT");
480  "DISPLACEMENT");
482  "DISPLACEMENT");
483
484  if (m_field.check_field("MESH_NODE_POSITIONS")) {
486  "DAMPER", "MESH_NODE_POSITIONS");
487  }
488  std::map<int, KelvinVoigtDamper::BlockMaterialData>::iterator bit =
489  damper.blockMaterialDataMap.begin();
490  for (; bit != damper.blockMaterialDataMap.end(); bit++) {
491  bit->second.lInear = linear;
492  int id = bit->first;
493  KelvinVoigtDamper::BlockMaterialData &material_data = bit->second;
494  damper.constitutiveEquationMap.insert(
496  material_data));
498  MBTET, "DAMPER");
499  }
500  CHKERR damper.setOperators(3);
501  }
502
503  MonitorPostProc post_proc(m_field, elastic.setOfBlocks,
504  elastic.getLoopFeEnergy(),
505  inertia.getLoopFeEnergy());
506
507  // elastic and mass element calculated in Kuu shell matrix problem. To
508  // calculate Mass element, velocity field is needed.
511  "DOT_DISPLACEMENT");
513  "DOT_VELOCITY");
514
515  // build field
516  CHKERR m_field.build_fields();
517  // CHKERR m_field.list_dofs_by_field_name("DISPLACEMENT");
518
519  // 10 node tets
520  if (!check_if_spatial_field_exist) {
521  Projection10NodeCoordsOnField ent_method_material(m_field,
522  "MESH_NODE_POSITIONS");
523  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
524  }
525
526  // build finite elements
527  CHKERR m_field.build_finite_elements();
530
531  // define problems
532  {
538  "FLUID_PRESSURE_FE");
540
541  ProblemsManager *prb_mng_ptr;
542  CHKERR m_field.getInterface(prb_mng_ptr);
543  if (is_partitioned) {
544  CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("Kuu", true);
545  CHKERR prb_mng_ptr->partitionFiniteElements("Kuu", true, 0,
546  pcomm->size());
547  } else {
548  CHKERR prb_mng_ptr->buildProblem("Kuu", true);
549  CHKERR prb_mng_ptr->partitionProblem("Kuu");
550  CHKERR prb_mng_ptr->partitionFiniteElements("Kuu");
551  }
552  CHKERR prb_mng_ptr->partitionGhostDofs("Kuu");
553  }
554
556  // set finite elements for problems
562  "FLUID_PRESSURE_FE");
564  "MASS_ELEMENT");
566  "VELOCITY_ELEMENT");
567  // set refinement level for problem
569
570  ProblemsManager *prb_mng_ptr;
571  CHKERR m_field.getInterface(prb_mng_ptr);
572  if (is_partitioned) {
573  CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("DYNAMICS", true);
574  CHKERR prb_mng_ptr->partitionFiniteElements("DYNAMICS", true, 0,
575  pcomm->size());
576  } else {
577  CHKERR prb_mng_ptr->buildProblem("DYNAMICS", true);
578  CHKERR prb_mng_ptr->partitionProblem("DYNAMICS");
579  CHKERR prb_mng_ptr->partitionFiniteElements("DYNAMICS");
580  }
581  CHKERR prb_mng_ptr->partitionGhostDofs("DYNAMICS");
582
583  Vec F;
584  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("DYNAMICS", COL,
585  &F);
586  Vec D;
587  CHKERR VecDuplicate(F, &D);
588
589  // create tS
590  TS ts;
591  CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
592  CHKERR TSSetType(ts, TSBEULER);
593
594  // shell matrix
595  ConvectiveMassElement::MatShellCtx *shellAij_ctx =
598  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("Kuu",
599  &shellAij_ctx->K);
600  CHKERR MatDuplicate(shellAij_ctx->K, MAT_DO_NOT_COPY_VALUES,
601  &shellAij_ctx->M);
602  CHKERR shellAij_ctx->iNit();
603  CHKERR m_field.getInterface<VecManager>()->vecScatterCreate(
604  D, "DYNAMICS", COL, shellAij_ctx->u, "Kuu", COL,
605  &shellAij_ctx->scatterU);
606  CHKERR m_field.getInterface<VecManager>()->vecScatterCreate(
607  D, "DYNAMICS", "VELOCITY", COL, shellAij_ctx->v, "Kuu", "DISPLACEMENT",
608  COL, &shellAij_ctx->scatterV);
609  Mat shell_Aij;
610  const Problem *problem_ptr;
611  CHKERR m_field.get_problem("DYNAMICS", &problem_ptr);
612  CHKERR MatCreateShell(
613  PETSC_COMM_WORLD, problem_ptr->getNbLocalDofsRow(),
614  problem_ptr->getNbLocalDofsCol(), problem_ptr->getNbDofsRow(),
615  problem_ptr->getNbDofsRow(), (void *)shellAij_ctx, &shell_Aij);
616  CHKERR MatShellSetOperation(shell_Aij, MATOP_MULT,
617  (void (*)(void))ConvectiveMassElement::MultOpA);
618  CHKERR MatShellSetOperation(
619  shell_Aij, MATOP_ZERO_ENTRIES,
620  (void (*)(void))ConvectiveMassElement::ZeroEntriesOp);
621  // blocked problem
622  ConvectiveMassElement::ShellMatrixElement shell_matrix_element(m_field);
623  DirichletDisplacementBc shell_dirichlet_bc(
624  m_field, "DISPLACEMENT", shellAij_ctx->barK, PETSC_NULL, PETSC_NULL);
625  DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT", PETSC_NULL,
626  D, F);
627  shell_matrix_element.problemName = "Kuu";
628  shell_matrix_element.shellMatCtx = shellAij_ctx;
629  shell_matrix_element.DirichletBcPtr = &shell_dirichlet_bc;
630  shell_matrix_element.loopK.push_back(
631  ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
632  "ELASTIC", &elastic.getLoopFeLhs()));
633  // damper
634  shell_matrix_element.loopK.push_back(
635  ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
636  "ELASTIC", &damper.feLhs));
637
639  CHKERR inertia.setShellMatrixMassOperators("VELOCITY", "DISPLACEMENT",
640  "MESH_NODE_POSITIONS", linear);
641  // element name "ELASTIC" is used, therefore M matrix is assembled as K
642  // matrix. This is added to M is shell matrix. M matrix is a derivative of
643  // inertia forces over spatial velocities
644  shell_matrix_element.loopM.push_back(
645  ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
646  "ELASTIC", &inertia.getLoopFeMassLhs()));
647  // this calculate derivatives of inertia forces over spatial positions and
648  // add this to shell K matrix
649  shell_matrix_element.loopAuxM.push_back(
650  ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
651  "ELASTIC", &inertia.getLoopFeMassAuxLhs()));
652
653  // Element to calculate shell matrix residual
654  ConvectiveMassElement::ShellResidualElement shell_matrix_residual(m_field);
655  shell_matrix_residual.shellMatCtx = shellAij_ctx;
656
657  // surface pressure
658  boost::ptr_map<std::string, NeumannForcesSurface> surface_forces;
659  {
660  string fe_name_str = "FORCE_FE";
661  surface_forces.insert(fe_name_str, new NeumannForcesSurface(m_field));
663  surface_forces.at(fe_name_str).getLoopFe(), false,
664  false);
666  NODESET | FORCESET, it)) {
667  CHKERR surface_forces.at(fe_name_str)
669  surface_forces.at(fe_name_str)
670  .methodsOp.push_back(new TimeForceScale());
671  }
672  }
673
674  boost::ptr_map<std::string, NeumannForcesSurface> surface_pressure;
675  {
676  string fe_name_str = "PRESSURE_FE";
677  surface_pressure.insert(fe_name_str, new NeumannForcesSurface(m_field));
679  surface_pressure.at(fe_name_str).getLoopFe(), false,
680  false);
682  m_field, SIDESET | PRESSURESET, it)) {
683  CHKERR surface_pressure.at(fe_name_str)
685  surface_pressure.at(fe_name_str)
686  .methodsOp.push_back(new TimeForceScale());
687  }
688  }
689
690  // edge forces
691  boost::ptr_map<std::string, EdgeForce> edge_forces;
692  {
693  string fe_name_str = "FORCE_FE";
694  edge_forces.insert(fe_name_str, new EdgeForce(m_field));
696  NODESET | FORCESET, it)) {
697  CHKERR edge_forces.at(fe_name_str)
699  edge_forces.at(fe_name_str).methodsOp.push_back(new TimeForceScale());
700  }
701  }
702
703  // nodal forces
704  boost::ptr_map<std::string, NodalForce> nodal_forces;
705  {
706  string fe_name_str = "FORCE_FE";
707  nodal_forces.insert(fe_name_str, new NodalForce(m_field));
709  NODESET | FORCESET, it)) {
710  CHKERR nodal_forces.at(fe_name_str)
712  nodal_forces.at(fe_name_str).methodsOp.push_back(new TimeForceScale());
713  }
714  }
715
716  MonitorRestart monitor_restart(m_field, ts);
717  ConvectiveMassElement::UpdateAndControl update_and_control(
718  m_field, ts, "VELOCITY", "DISPLACEMENT");
719
720  // TS
721  TsCtx ts_ctx(m_field, "DYNAMICS");
722
723  // right hand side
724  // preprocess
725  ts_ctx.getPreProcessIFunction().push_back(&update_and_control);
726  ts_ctx.getPreProcessIFunction().push_back(&my_dirichlet_bc);
727
728  // fe looops
729  auto &loops_to_do_Rhs = ts_ctx.getLoopsIFunction();
730
731  auto add_static_rhs = [&](auto &loops_to_do_Rhs) {
733  loops_to_do_Rhs.push_back(
734  PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeRhs()));
735  for (auto fit = surface_forces.begin(); fit != surface_forces.end();
736  fit++) {
737  loops_to_do_Rhs.push_back(
738  PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
739  }
740  for (auto fit = surface_pressure.begin(); fit != surface_pressure.end();
741  fit++) {
742  loops_to_do_Rhs.push_back(
743  PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
744  }
745  for (auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
746  loops_to_do_Rhs.push_back(
747  PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
748  }
749  for (auto fit = nodal_forces.begin(); fit != nodal_forces.end(); fit++) {
750  loops_to_do_Rhs.push_back(
751  PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
752  }
753  loops_to_do_Rhs.push_back(PairNameFEMethodPtr(
754  "FLUID_PRESSURE_FE", &fluid_pressure_fe.getLoopFe()));
756  };
757
759
760  loops_to_do_Rhs.push_back(PairNameFEMethodPtr("DAMPER", &damper.feRhs));
761  loops_to_do_Rhs.push_back(
762  PairNameFEMethodPtr("MASS_ELEMENT", &inertia.getLoopFeMassRhs()));
763
764  // preporcess
765  // calculate residual for velocities
766  ts_ctx.getPreProcessIFunction().push_back(&shell_matrix_residual);
767  // postprocess
768  ts_ctx.getPostProcessIFunction().push_back(&my_dirichlet_bc);
769
770  // left hand side
771  // preprocess
772  ts_ctx.getPreProcessIJacobian().push_back(&update_and_control);
773  ts_ctx.getPreProcessIJacobian().push_back(&shell_matrix_element);
774  ts_ctx.getPostProcessIJacobian().push_back(&update_and_control);
775  // monitor
776  TsCtx::FEMethodsSequence &loopsMonitor =
778  loopsMonitor.push_back(
779  TsCtx::PairNameFEMethodPtr("MASS_ELEMENT", &post_proc));
780  loopsMonitor.push_back(
781  TsCtx::PairNameFEMethodPtr("MASS_ELEMENT", &monitor_restart));
782
783  CHKERR TSSetIFunction(ts, F, TsSetIFunction, &ts_ctx);
784  CHKERR TSSetIJacobian(ts, shell_Aij, shell_Aij, TsSetIJacobian, &ts_ctx);
785
786  CHKERR TSMonitorSet(ts, TsMonitorSet, &ts_ctx, PETSC_NULL);
787
788  double ftime = 1;
789  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
790  CHKERR TSSetSolution(ts, D);
791  CHKERR TSSetFromOptions(ts);
792  // shell matrix pre-conditioner
793  SNES snes;
794  CHKERR TSGetSNES(ts, &snes);
795  // CHKERR SNESSetFromOptions(snes);
796  KSP ksp;
797  CHKERR SNESGetKSP(snes, &ksp);
798  CHKERR KSPSetFromOptions(ksp);
799  PC pc;
800  CHKERR KSPGetPC(ksp, &pc);
801  CHKERR PCSetType(pc, PCSHELL);
802  ConvectiveMassElement::PCShellCtx pc_shell_ctx(shell_Aij);
803  CHKERR PCShellSetContext(pc, (void *)&pc_shell_ctx);
804  CHKERR PCShellSetApply(pc, ConvectiveMassElement::PCShellApplyOp);
805  CHKERR PCShellSetSetUp(pc, ConvectiveMassElement::PCShellSetUpOp);
806  CHKERR PCShellSetDestroy(pc, ConvectiveMassElement::PCShellDestroy);
807
808  CHKERR VecZeroEntries(D);
809  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
810  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
811  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
812  "DYNAMICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
813
814  // Solve problem at time Zero
815  if (is_solve_at_time_zero) {
816
817  Mat Aij = shellAij_ctx->K;
818  Vec F;
819  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("Kuu", COL, &F);
820  Vec D;
821  CHKERR VecDuplicate(F, &D);
822
823  // Set vector for Kuu problem from the mesh data
824  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
825  "Kuu", COL, D, INSERT_VALUES, SCATTER_FORWARD);
826  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
827  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
828
829  SnesCtx snes_ctx(m_field, "Kuu");
830
831  SNES snes;
832  CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
833  CHKERR SNESSetApplicationContext(snes, &snes_ctx);
834  CHKERR SNESSetFunction(snes, F, SnesRhs, &snes_ctx);
835  CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, &snes_ctx);
836  CHKERR SNESSetFromOptions(snes);
837
838  DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT",
839  PETSC_NULL, D, F);
840
841  SnesCtx::FEMethodsSequence &loops_to_do_Rhs =
842  snes_ctx.getComputeRhs();
843  snes_ctx.getPreProcComputeRhs().push_back(&my_dirichlet_bc);
844  fluid_pressure_fe.getLoopFe().ts_t = 0;
846  snes_ctx.getPostProcComputeRhs().push_back(&my_dirichlet_bc);
847
848  SnesCtx::FEMethodsSequence &loops_to_do_Mat =
849  snes_ctx.getSetOperators();
850  snes_ctx.getPreProcSetOperators().push_back(&my_dirichlet_bc);
851  loops_to_do_Mat.push_back(
852  SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeLhs()));
853  snes_ctx.getPostProcSetOperators().push_back(&my_dirichlet_bc);
854
855  CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0, "VELOCITY");
856  CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0,
857  "DOT_DISPLACEMENT");
858  CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0, "DOT_VELOCITY");
859
860  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
861  "Kuu", COL, D, INSERT_VALUES, SCATTER_FORWARD);
862
863  CHKERR SNESSolve(snes, PETSC_NULL, D);
864  int its;
865  CHKERR SNESGetIterationNumber(snes, &its);
866  MOFEM_LOG_C("DYNAMIC", Sev::inform, "number of Newton iterations = %d\n",
867  its);
868
869  // Set data on the mesh
870  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
871  "Kuu", COL, D, INSERT_VALUES, SCATTER_REVERSE);
872
873  CHKERR VecDestroy(&F);
874  CHKERR VecDestroy(&D);
875  CHKERR SNESDestroy(&snes);
876  }
877
878  if (is_solve_at_time_zero) {
879  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
880  "DYNAMICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
881  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
882  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
883  CHKERR TSSetSolution(ts, D);
884  }
885
886 #if PETSC_VERSION_GE(3, 7, 0)
887  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
888 #endif
889  CHKERR TSSolve(ts, D);
890  CHKERR TSGetTime(ts, &ftime);
891
892  PetscInt steps, snesfails, rejects, nonlinits, linits;
893  CHKERR TSGetTimeStepNumber(ts, &steps);
894  CHKERR TSGetSNESFailures(ts, &snesfails);
895  CHKERR TSGetStepRejections(ts, &rejects);
896  CHKERR TSGetSNESIterations(ts, &nonlinits);
897  CHKERR TSGetKSPIterations(ts, &linits);
898  MOFEM_LOG_C("DYNAMIC", Sev::inform,
899  "steps %d (%d rejected, %D SNES fails), ftime %g, nonlinits "
900  "%d, linits %D\n",
901  steps, rejects, snesfails, ftime, nonlinits, linits);
902  CHKERR TSDestroy(&ts);
903
904  CHKERR VecDestroy(&F);
905  CHKERR VecDestroy(&D);
906  CHKERR MatDestroy(&shellAij_ctx->K);
907  CHKERR MatDestroy(&shellAij_ctx->M);
908  CHKERR VecScatterDestroy(&shellAij_ctx->scatterU);
909  CHKERR VecScatterDestroy(&shellAij_ctx->scatterV);
910  CHKERR MatDestroy(&shell_Aij);
911  delete shellAij_ctx;
912  }
913  CATCH_ERRORS;
914
916
917  return 0;
918 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::SnesCtx::getPostProcComputeRhs
BasicMethodsSequence & getPostProcComputeRhs()
Definition: SnesCtx.hpp:87
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
Definition: HODataOperators.hpp:789
SIDESET
@ SIDESET
Definition: definitions.h:160
MoFEM::CoreInterface::loop_finite_elements
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEMErrorCode addNeumannFluidPressureBCElements(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Definition: FluidPressure.cpp:67
static MoFEMErrorCode addNeumannBCElements(MoFEM::Interface &m_field, const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS", Range *intersect_ptr=NULL)
Declare finite element.
Definition: SurfacePressure.cpp:1974
MonitorPostProc::mField
MoFEM::Interface & mField
Definition: nonlinear_dynamics.cpp:32
MoFEM::CoreInterface::loop_dofs
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
MoFEM::ProblemsManager::buildProblem
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
Definition: ProblemsManager.cpp:87
KelvinVoigtDamper::CommonData::spatialPositionNameDot
string spatialPositionNameDot
Definition: KelvinVoigtDamper.hpp:172
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
EdgeForce
Force on edges and lines.
Definition: EdgeForce.hpp:13
MonitorPostProc::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: nonlinear_dynamics.cpp:133
ConvectiveMassElement::MatShellCtx::iNit
MoFEMErrorCode iNit()
Definition: ConvectiveMassElement.cpp:2352
ConvectiveMassElement::setShellMatrixMassOperators
MoFEMErrorCode setShellMatrixMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool linear=false)
Definition: ConvectiveMassElement.cpp:2243
EntityHandle
MonitorPostProc::feElasticEnergy
NonlinearElasticElement::MyVolumeFE & feElasticEnergy
calculate elastic energy
Definition: nonlinear_dynamics.cpp:36
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:92
PostProcTemplateOnRefineMesh::postProcMesh
moab::Interface & postProcMesh
Definition: PostProcOnRefMesh.hpp:122
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
Definition: HODataOperators.hpp:764
PRESSURESET
@ PRESSURESET
Definition: definitions.h:165
MoFEM::TsSetIFunction
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
Definition: TsCtx.cpp:56
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
NOBASE
@ NOBASE
Definition: definitions.h:59
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
NonlinearElasticElement
structure grouping operators and data used for calculation of nonlinear elastic element
Definition: HookeElement.hpp:27
MoFEM::TsCtx::getLoopsMonitor
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition: TsCtx.hpp:98
KelvinVoigtDamper
Implementation of Kelvin Voigt Damper.
Definition: KelvinVoigtDamper.hpp:20
MoFEM::TsCtx::getLoopsIFunction
FEMethodsSequence & getLoopsIFunction()
Get the loops to do IFunction object.
Definition: TsCtx.hpp:59
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
E
_IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
Definition: MeshsetsManager.hpp:49
MoFEM::SnesCtx::getPreProcSetOperators
BasicMethodsSequence & getPreProcSetOperators()
Definition: SnesCtx.hpp:92
MonitorRestart::mField
MoFEM::Interface & mField
Definition: nonlinear_dynamics.cpp:143
MoFEM::CoreInterface::get_problem
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
NodalForce
Force applied to nodes.
Definition: NodalForce.hpp:13
ConvectiveMassElement::MatShellCtx::scatterU
VecScatter scatterU
Definition: ConvectiveMassElement.hpp:505
MoFEM::TsCtx::getPostProcessIJacobian
BasicMethodsSequence & getPostProcessIJacobian()
Get the postProcess to do IJacobian object.
Definition: TsCtx.hpp:128
BasicFiniteElements.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
KelvinVoigtDamper::feLhs
DamperFE feLhs
Definition: KelvinVoigtDamper.hpp:247
ElasticMaterials
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
Definition: ElasticMaterials.hpp:24
MoFEM::TsCtx::FEMethodsSequence
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: TsCtx.hpp:26
MoFEM::SnesCtx::FEMethodsSequence
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: SnesCtx.hpp:27
MonitorPostProc::iNit
bool iNit
Definition: nonlinear_dynamics.cpp:40
MoFEM::LogManager::createSink
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:298
ts_ctx
MoFEM::TsCtx * ts_ctx
Definition: level_set.cpp:1932
DirichletDisplacementBc
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:57
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::SnesCtx::getPostProcSetOperators
BasicMethodsSequence & getPostProcSetOperators()
Definition: SnesCtx.hpp:101
MoFEMErrorCode addVelocityElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
Definition: ConvectiveMassElement.cpp:1892
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::TsMonitorSet
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:259
ConvectiveMassElement
structure grouping operators and data used for calculation of mass (convective) element
Definition: ConvectiveMassElement.hpp:28
NODESET
@ NODESET
Definition: definitions.h:159
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
MoFEM::Exceptions::rval
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:74
ConvectiveMassElement::MatShellCtx::u
Vec u
Definition: ConvectiveMassElement.hpp:513
MonitorRestart::MonitorRestart
MonitorRestart(MoFEM::Interface &m_field, TS ts)
Definition: nonlinear_dynamics.cpp:146
SurfacePressureComplexForLazy.hpp
ConvectiveMassElement::MyVolumeFE
definition of volume element
Definition: ConvectiveMassElement.hpp:31
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEMErrorCode addConvectiveMassElement(string element_name, string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool ale=false, BitRefLevel bit=BitRefLevel())
Definition: ConvectiveMassElement.cpp:1835
ConvectiveMassElement::PCShellSetUpOp
static MoFEMErrorCode PCShellSetUpOp(PC pc)
Definition: ConvectiveMassElement.hpp:618
MonitorRestart::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: nonlinear_dynamics.cpp:217
PostProcStress
Definition: PostProcStresses.hpp:17
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
KelvinVoigtDamper::setOperators
MoFEMErrorCode setOperators(const int tag)
Definition: KelvinVoigtDamper.hpp:812
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
MonitorPostProc::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: nonlinear_dynamics.cpp:128
ConvectiveMassElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: ConvectiveMassElement.hpp:122
FORCESET
@ FORCESET
Definition: definitions.h:164
TimeForceScale
Force scale operator for reading two columns.
Definition: TimeForceScale.hpp:18
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
KelvinVoigtDamper::blockMaterialDataMap
std::map< int, BlockMaterialData > blockMaterialDataMap
Definition: KelvinVoigtDamper.hpp:37
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
MoFEM::TsCtx::getPreProcessIJacobian
BasicMethodsSequence & getPreProcessIJacobian()
Get the preProcess to do IJacobian object.
Definition: TsCtx.hpp:121
ConvectiveMassElement::PCShellDestroy
static MoFEMErrorCode PCShellDestroy(PC pc)
Definition: ConvectiveMassElement.hpp:633
MoFEM::SnesMethod::CTX_SNESNONE
@ CTX_SNESNONE
Definition: LoopMethods.hpp:107
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
MoFEM::Problem::getNbDofsRow
DofIdx getNbDofsRow() const
Definition: ProblemsMultiIndices.hpp:376
PostProcTemplateVolumeOnRefinedMesh::commonData
CommonData commonData
Definition: PostProcOnRefMesh.hpp:287
TimeForceScale.hpp
FluidPressure
Fluid pressure forces.
Definition: FluidPressure.hpp:20
MonitorRestart::step
int * step
Definition: nonlinear_dynamics.cpp:142
MoFEM::PairNameFEMethodPtr
Definition: AuxPETSc.hpp:12
COL
@ COL
Definition: definitions.h:136
ConvectiveMassElement::ShellResidualElement::shellMatCtx
MatShellCtx * shellMatCtx
pointer to shell matrix
Definition: ConvectiveMassElement.hpp:720
MoFEM::LogManager::getStrmWorld
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:344
PostProcTemplateOnRefineMesh::writeFile
MoFEMErrorCode writeFile(const std::string file_name, const char *file_type="MOAB", const char *file_options="PARALLEL=WRITE_PART")
wrote results in (MOAB) format, use "file_name.h5m"
Definition: PostProcOnRefMesh.hpp:253
MoFEM::SnesCtx::getComputeRhs
FEMethodsSequence & getComputeRhs()
Definition: SnesCtx.hpp:69
MoFEM::TsCtx
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:17
MoFEM::TsCtx::getPostProcessIFunction
BasicMethodsSequence & getPostProcessIFunction()
Get the postProcess to do IFunction object.
Definition: TsCtx.hpp:112
MoFEM::TSMethod::CTX_TSNONE
@ CTX_TSNONE
Definition: LoopMethods.hpp:145
ConvectiveMassElement::MatShellCtx::M
Mat M
Definition: ConvectiveMassElement.hpp:504
MonitorPostProc::MonitorPostProc
MonitorPostProc(MoFEM::Interface &m_field, std::map< int, NonlinearElasticElement::BlockData > &set_of_blocks, NonlinearElasticElement::MyVolumeFE &fe_elastic_energy, ConvectiveMassElement::MyVolumeFE &fe_kinetic_energy)
Definition: nonlinear_dynamics.cpp:45
main
int main(int argc, char *argv[])
Definition: nonlinear_dynamics.cpp:226
MonitorPostProc::setOfBlocks
std::map< int, NonlinearElasticElement::BlockData > & setOfBlocks
Definition: nonlinear_dynamics.cpp:34
ConvectiveMassElement::MatShellCtx::v
Vec v
Definition: ConvectiveMassElement.hpp:513
MoFEM::TsSetIJacobian
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobian in TS solver.
Definition: TsCtx.cpp:165
AINSWORTH_LOBATTO_BASE
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:62
PostProcTemplateVolumeOnRefinedMesh::generateReferenceElementMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
Definition: PostProcOnRefMesh.hpp:301
Definition: PostProcOnRefMesh.hpp:195
MoFEMErrorCode addElement(const std::string element_name, const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false)
Definition: NonLinearElasticElement.cpp:1120
MoFEM::CoreInterface::check_field
virtual bool check_field(const std::string &name) const =0
check if field is in database
ConvectiveMassElement::ZeroEntriesOp
static MoFEMErrorCode ZeroEntriesOp(Mat A)
Definition: ConvectiveMassElement.hpp:589
ConvectiveMassElement::MatShellCtx::K
Mat K
Definition: ConvectiveMassElement.hpp:504
ConvectiveMassElement::PCShellApplyOp
static MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x)
apply pre-conditioner for shell matrix
Definition: ConvectiveMassElement.hpp:671
MOFEM_LOG_TAG
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
MonitorPostProc::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: nonlinear_dynamics.cpp:80
MoFEM::Problem::getNbLocalDofsRow
DofIdx getNbLocalDofsRow() const
Definition: ProblemsMultiIndices.hpp:378
ElasticMaterials.hpp
Elastic materials.
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::ProblemsManager::partitionFiniteElements
MoFEMErrorCode partitionFiniteElements(const std::string name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=VERBOSE)
partition finite elements
Definition: ProblemsManager.cpp:2167
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
MoFEM::VecManager
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:23
NonlinearElasticElement::setOperators
MoFEMErrorCode setOperators(const std::string spatial_position_field_name, const std::string material_position_field_name="MESH_NODE_POSITIONS", const bool ale=false, const bool field_disp=false)
Set operators to calculate left hand tangent matrix and right hand residual.
Definition: NonLinearElasticElement.cpp:1153
AINSWORTH_BERNSTEIN_BEZIER_BASE
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
mesh_file_name
char mesh_file_name[255]
Definition: mesh_smoothing.cpp:23
MoFEM::TsCtx::getPreProcessIFunction
BasicMethodsSequence & getPreProcessIFunction()
Get the preProcess to do IFunction object.
Definition: TsCtx.hpp:105
MoFEM::MatrixManager
Matrix manager is used to build and partition problems.
Definition: MatrixManager.hpp:21
MoFEM::SnesCtx::getSetOperators
FEMethodsSequence & getSetOperators()
Definition: SnesCtx.hpp:64
ConvectiveMassElement::methodsOp
boost::ptr_vector< MethodForForceScaling > methodsOp
Definition: ConvectiveMassElement.hpp:150
MonitorRestart::time
double * time
Definition: nonlinear_dynamics.cpp:141
static MoFEMErrorCode addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: EdgeForce.hpp:62
NonlinearElasticElement::MyVolumeFE
definition of volume element
Definition: NonLinearElasticElement.hpp:31
MoFEM::SnesRhs
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
Definition: SnesCtx.cpp:27
MoFEM::SnesCtx::getPreProcComputeRhs
BasicMethodsSequence & getPreProcComputeRhs()
Definition: SnesCtx.hpp:78
MonitorPostProc::step
int * step
Definition: nonlinear_dynamics.cpp:43
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
FluidPressure::getLoopFe
MyTriangleFE & getLoopFe()
Definition: FluidPressure.hpp:35
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::TSMethod::ts_t
PetscReal ts_t
time
Definition: LoopMethods.hpp:162
PostProcVolumeOnRefinedMesh
Post processing.
Definition: PostProcOnRefMesh.hpp:955
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
KelvinVoigtDamper::feRhs
DamperFE feRhs
Definition: KelvinVoigtDamper.hpp:247
FluidPressure::setNeumannFluidPressureFiniteElementOperators
MoFEMErrorCode setNeumannFluidPressureFiniteElementOperators(string field_name, Vec F, bool allow_negative_pressure=true, bool ho_geometry=false)
Definition: FluidPressure.cpp:131
PostProcTemplateOnRefineMesh::mapGaussPts
std::vector< EntityHandle > mapGaussPts
Definition: PostProcOnRefMesh.hpp:125
ConvectiveMassElement::UpdateAndControl
Set fields DOT_.
Definition: ConvectiveMassElement.hpp:434
MonitorRestart
Definition: nonlinear_dynamics.cpp:139
Definition: ConvectiveMassElement.hpp:92
MonitorRestart::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: nonlinear_dynamics.cpp:212
NonlinearElasticElement::setOfBlocks
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
Definition: NonLinearElasticElement.hpp:100
MonitorPostProc::feKineticEnergy
ConvectiveMassElement::MyVolumeFE & feKineticEnergy
calculate elastic energy
Definition: nonlinear_dynamics.cpp:38
ConvectiveMassElement::MultOpA
static MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f)
Mult operator for shell matrix.
Definition: ConvectiveMassElement.hpp:546
KelvinVoigtDamper::CommonData
Common data for nonlinear_elastic_elem model.
Definition: KelvinVoigtDamper.hpp:169
MoFEM::Exceptions::ierr
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Exceptions.hpp:76
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::ForcesAndSourcesCore::getOpPtrVector
boost::ptr_deque< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
Definition: ForcesAndSourcesCore.hpp:83
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
MoFEM::ProblemsManager::buildProblemOnDistributedMesh
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)
Definition: ProblemsManager.cpp:290
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
MoFEM::PetscOptionsGetEList
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
Definition: DeprecatedPetsc.hpp:203
MoFEM::CoreInterface::build_fields
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
FieldApproximationBase
FieldApproximationBase
approximation base
Definition: definitions.h:58
MoFEM::Problem::getNbLocalDofsCol
DofIdx getNbLocalDofsCol() const
Definition: ProblemsMultiIndices.hpp:379
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
help
static char help[]
Definition: nonlinear_dynamics.cpp:28
MoFEM::BitRefManager
Managing BitRefLevels.
Definition: BitRefManager.hpp:21
virtual MoFEMErrorCode modify_problem_add_finite_element(const std::string name_problem, const std::string &fe_name)=0
add finite element to problem, this add entities assigned to finite element to a particular problem
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:40
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
ConvectiveMassElement::MatShellCtx::scatterV
VecScatter scatterV
Definition: ConvectiveMassElement.hpp:505
MoFEM::Problem
Definition: ProblemsMultiIndices.hpp:54
KelvinVoigtDamper::CommonData::spatialPositionName
string spatialPositionName
Definition: KelvinVoigtDamper.hpp:171
MonitorRestart::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: nonlinear_dynamics.cpp:195
MoFEM::CoreInterface::set_field_order
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEM::LogManager::setLog
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:389
NeumannForcesSurface
Finite element and operators to apply force/pressures applied to surfaces.
Definition: SurfacePressure.hpp:14
MonitorPostProc
Definition: nonlinear_dynamics.cpp:30
KelvinVoigtDamper::constitutiveEquationMap
ConstitutiveEquationMap constitutiveEquationMap
Definition: KelvinVoigtDamper.hpp:164
MoFEM::SnesMat
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
Definition: SnesCtx.cpp:139
MoFEM::ProblemsManager::partitionGhostDofs
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
Definition: ProblemsManager.cpp:2339
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
MonitorRestart::pRT
int pRT
Definition: nonlinear_dynamics.cpp:144
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
TimeAccelerogram
Definition: TimeForceScale.hpp:133
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
virtual MoFEMErrorCode add_field(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MoFEM::ProblemsManager::partitionProblem
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
Definition: ProblemsManager.cpp:1683
MonitorPostProc::pRT
int pRT
Definition: nonlinear_dynamics.cpp:42
F
@ F
Definition: free_surface.cpp:394
KelvinVoigtDamper::commonData
CommonData commonData
Definition: KelvinVoigtDamper.hpp:198
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182