v0.13.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  * This file is part of MoFEM.
22  * MoFEM is free software: you can redistribute it and/or modify it under
23  * the terms of the GNU Lesser General Public License as published by the
24  * Free Software Foundation, either version 3 of the License, or (at your
25  * option) any later version.
26  *
27  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
28  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
29  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
30  * License for more details.
31  *
32  * You should have received a copy of the GNU Lesser General Public
33  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
34 
35 #include <BasicFiniteElements.hpp>
36 using namespace MoFEM;
37 
38 #include <ElasticMaterials.hpp>
40 
41 static char help[] = "...\n\n";
42 
43 struct MonitorPostProc : public FEMethod {
44 
47  std::map<int, NonlinearElasticElement::BlockData> &setOfBlocks;
49  &feElasticEnergy; ///< calculate elastic energy
51  &feKineticEnergy; ///< calculate elastic energy
52 
53  bool iNit;
54 
55  int pRT;
56  int *step;
57 
59  MoFEM::Interface &m_field,
60  std::map<int, NonlinearElasticElement::BlockData> &set_of_blocks,
61  NonlinearElasticElement::MyVolumeFE &fe_elastic_energy,
62  ConvectiveMassElement::MyVolumeFE &fe_kinetic_energy)
63  : FEMethod(), mField(m_field), postProc(m_field),
64  setOfBlocks(set_of_blocks), feElasticEnergy(fe_elastic_energy),
65  feKineticEnergy(fe_kinetic_energy), iNit(false) {
66 
67  double def_t_val = 0;
68  const EntityHandle root_meshset = mField.get_moab().get_root_set();
69 
70  Tag th_step;
71  rval = m_field.get_moab().tag_get_handle(
72  "_TsStep_", 1, MB_TYPE_INTEGER, th_step,
73  MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_t_val);
74  if (rval == MB_ALREADY_ALLOCATED) {
75  MOAB_THROW(m_field.get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
76  (const void **)&step));
77  } else {
78  MOAB_THROW(m_field.get_moab().tag_set_data(th_step, &root_meshset, 1,
79  &def_t_val));
80  MOAB_THROW(m_field.get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
81  (const void **)&step));
82  }
83 
84  PetscBool flg = PETSC_TRUE;
85  CHK_THROW_MESSAGE(PetscOptionsGetInt(PETSC_NULL, PETSC_NULL,
86  "-my_output_prt", &pRT, &flg),
87  "Can not get option");
88  if (flg != PETSC_TRUE) {
89  pRT = 10;
90  }
91  }
92 
95 
96  if (!iNit) {
98  if(mField.check_field("MESH_NODE_POSITIONS"))
99  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", postProc, true, false, false,
100  false);
101  CHKERR postProc.addFieldValuesPostProc("DISPLACEMENT");
102  CHKERR postProc.addFieldValuesPostProc("VELOCITY");
103  CHKERR postProc.addFieldValuesPostProc("MESH_NODE_POSITIONS");
104  CHKERR postProc.addFieldValuesGradientPostProc("DISPLACEMENT");
105 
106  std::map<int, NonlinearElasticElement::BlockData>::iterator sit =
107  setOfBlocks.begin();
108  for (; sit != setOfBlocks.end(); sit++) {
109  postProc.getOpPtrVector().push_back(new PostProcStress(
110  postProc.postProcMesh, postProc.mapGaussPts, "DISPLACEMENT",
111  sit->second, postProc.commonData, true));
112  }
113 
114  iNit = true;
115  }
116 
117  if ((*step) % pRT == 0) {
118  CHKERR mField.loop_finite_elements("DYNAMICS", "MASS_ELEMENT", postProc);
119  std::ostringstream sss;
120  sss << "out_values_" << (*step) << ".h5m";
121  CHKERR postProc.writeFile(sss.str().c_str());
122  }
123 
124  feElasticEnergy.ts_ctx = TSMethod::CTX_TSNONE;
125  feElasticEnergy.snes_ctx = SnesMethod::CTX_SNESNONE;
126  CHKERR mField.loop_finite_elements("DYNAMICS", "ELASTIC", feElasticEnergy);
127 
128  feKineticEnergy.ts_ctx = TSMethod::CTX_TSNONE;
129  CHKERR mField.loop_finite_elements("DYNAMICS", "MASS_ELEMENT",
130  feKineticEnergy);
131  double E = feElasticEnergy.eNergy;
132  double T = feKineticEnergy.eNergy;
133  MOFEM_LOG_C(
134  "DYNAMIC", Sev::inform,
135  "%d Time %3.2e Elastic energy %3.2e Kinetic Energy %3.2e Total %3.2e\n",
136  ts_step, ts_t, E, T, E + T);
137 
139  }
140 
144  }
145 
149  }
150 };
151 
152 struct MonitorRestart : public FEMethod {
153 
154  double *time;
155  int *step;
157  int pRT;
158 
159  MonitorRestart(MoFEM::Interface &m_field, TS ts) : mField(m_field) {
160  double def_t_val = 0;
161 
162  const EntityHandle root_meshset = mField.get_moab().get_root_set();
163 
164  Tag th_time;
165  rval = m_field.get_moab().tag_get_handle(
166  "_TsTime_", 1, MB_TYPE_DOUBLE, th_time,
167  MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_t_val);
168  if (rval == MB_ALREADY_ALLOCATED) {
169  rval = m_field.get_moab().tag_get_by_ptr(th_time, &root_meshset, 1,
170  (const void **)&time);
171  MOAB_THROW(rval);
172  ierr = TSSetTime(ts, *time);
173  CHKERRABORT(PETSC_COMM_WORLD, ierr);
174  } else {
175  rval = m_field.get_moab().tag_set_data(th_time, &root_meshset, 1,
176  &def_t_val);
177  MOAB_THROW(rval);
178  rval = m_field.get_moab().tag_get_by_ptr(th_time, &root_meshset, 1,
179  (const void **)&time);
180  MOAB_THROW(rval);
181  }
182  Tag th_step;
183  rval = m_field.get_moab().tag_get_handle(
184  "_TsStep_", 1, MB_TYPE_INTEGER, th_step,
185  MB_TAG_CREAT | MB_TAG_EXCL | MB_TAG_MESH, &def_t_val);
186  if (rval == MB_ALREADY_ALLOCATED) {
187  rval = m_field.get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
188  (const void **)&step);
189  MOAB_THROW(rval);
190  } else {
191  rval = m_field.get_moab().tag_set_data(th_step, &root_meshset, 1,
192  &def_t_val);
193  MOAB_THROW(rval);
194  rval = m_field.get_moab().tag_get_by_ptr(th_step, &root_meshset, 1,
195  (const void **)&step);
196  MOAB_THROW(rval);
197  }
198 
199  PetscBool flg = PETSC_TRUE;
200  ierr = PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_output_prt", &pRT,
201  &flg);
202  CHKERRABORT(PETSC_COMM_WORLD, ierr);
203  if (flg != PETSC_TRUE) {
204  pRT = 10;
205  }
206  }
207 
210 
211  (*time) = ts_t;
212  // if(pRT>0) {
213  // if((*step)%pRT==0) {
214  // std::ostringstream ss;
215  // ss << "restart_" << (*step) << ".h5m";
216  // CHKERR
217  // mField.get_moab().write_file(ss.str().c_str()/*,"MOAB","PARALLEL=WRITE_PART"*/);
218  //
219  // }
220  // }
221  (*step)++;
223  }
224 
228  }
229 
233  }
234 };
235 
236 // See file users_modules/elasticity/TimeForceScale.hpp
237 #include <TimeForceScale.hpp>
238 
239 int main(int argc, char *argv[]) {
240 
241  const string default_options = "-ksp_type fgmres \n"
242  "-pc_type lu \n"
243  "-pc_factor_mat_solver_type mumps \n"
244  "-mat_mumps_icntl_20 0 \n"
245  "-ksp_atol 1e-10 \n"
246  "-ksp_rtol 1e-10 \n"
247  "-snes_monitor \n"
248  "-snes_type newtonls \n"
249  "-snes_linesearch_type basic \n"
250  "-snes_max_it 100 \n"
251  "-snes_atol 1e-7 \n"
252  "-snes_rtol 1e-7 \n"
253  "-ts_monitor \n"
254  "-ts_type alpha \n";
255 
256  string param_file = "param_file.petsc";
257  if (!static_cast<bool>(ifstream(param_file))) {
258  std::ofstream file(param_file.c_str(), std::ios::ate);
259  if (file.is_open()) {
260  file << default_options;
261  file.close();
262  }
263  }
264 
265  MoFEM::Core::Initialize(&argc, &argv, param_file.c_str(), help);
266 
267  // Add logging channel for example
268  auto core_log = logging::core::get();
269  core_log->add_sink(
271  LogManager::setLog("DYNAMIC");
272  MOFEM_LOG_TAG("DYNAMIC", "dynamic");
273 
274  try {
275 
276  moab::Core mb_instance;
277  moab::Interface &moab = mb_instance;
278 
279  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
280  auto moab_comm_wrap =
281  boost::make_shared<WrapMPIComm>(PETSC_COMM_WORLD, false);
282  if (pcomm == NULL)
283  pcomm = new ParallelComm(&moab, moab_comm_wrap->get_comm());
284 
286  char mesh_file_name[255];
287  PetscBool is_partitioned = PETSC_FALSE;
288  PetscBool linear = PETSC_TRUE;
289  PetscInt disp_order = 1;
290  PetscInt vel_order = 1;
291  PetscBool is_solve_at_time_zero = PETSC_FALSE;
292 
293  auto read_command_line_parameters = [&]() {
295  PetscBool flg = PETSC_TRUE;
296  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, "-my_file",
297  mesh_file_name, 255, &flg);
298  if (flg != PETSC_TRUE)
299  SETERRQ(PETSC_COMM_SELF, 1, "Error -my_file (mesh file needed)");
300 
301  // use this if your mesh is partitioned and you run code on parts,
302  // you can solve very big problems
303  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-my_is_partitioned",
304  &is_partitioned, &flg);
305 
306  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL, "-is_linear", &linear,
307  PETSC_NULL);
308 
309  enum bases { LEGENDRE, LOBATTO, BERNSTEIN_BEZIER, LASBASETOP };
310  const char *list_bases[] = {"legendre", "lobatto", "bernstein_bezier"};
311  PetscInt choice_base_value = BERNSTEIN_BEZIER;
312  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
313  LASBASETOP, &choice_base_value, PETSC_NULL);
314  if (choice_base_value == LEGENDRE)
316  else if (choice_base_value == LOBATTO)
317  base = AINSWORTH_LOBATTO_BASE;
318  else if (choice_base_value == BERNSTEIN_BEZIER)
320 
321  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_disp_order",
322  &disp_order, &flg);
323  if (flg != PETSC_TRUE)
324  disp_order = 1;
325 
326  CHKERR PetscOptionsGetInt(PETSC_NULL, PETSC_NULL, "-my_vel_order",
327  &vel_order, &flg);
328  if (flg != PETSC_TRUE)
329  vel_order = disp_order;
330 
331  CHKERR PetscOptionsGetBool(PETSC_NULL, PETSC_NULL,
332  "-my_solve_at_time_zero",
333  &is_solve_at_time_zero, &flg);
334 
336  };
337 
338  auto read_mesh = [&]() {
340  if (is_partitioned == PETSC_TRUE) {
341  // Read mesh to MOAB
342  const char *option;
343  option = "PARALLEL=BCAST_DELETE;"
344  "PARALLEL_RESOLVE_SHARED_ENTS;"
345  "PARTITION=PARALLEL_PARTITION;";
346  CHKERR moab.load_file(mesh_file_name, 0, option);
347  } else {
348  const char *option;
349  option = "";
350  CHKERR moab.load_file(mesh_file_name, 0, option);
351  }
353  };
354 
355  CHKERR read_command_line_parameters();
356  CHKERR read_mesh();
357 
358  MoFEM::Core core(moab);
359  MoFEM::Interface &m_field = core;
360 
361  // ref meshset ref level 0
362  BitRefLevel bit_level0;
363  bit_level0.set(0);
364  EntityHandle meshset_level0;
365  CHKERR moab.create_meshset(MESHSET_SET, meshset_level0);
366  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevelByDim(
367  0, 3, bit_level0);
368  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
369  bit_level0, BitRefLevel().set(), meshset_level0);
370 
371  // Fields
372  CHKERR m_field.add_field("MESH_NODE_POSITIONS", H1, AINSWORTH_LEGENDRE_BASE,
373  3, MB_TAG_SPARSE, MF_ZERO);
374  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "MESH_NODE_POSITIONS");
375  CHKERR m_field.set_field_order(0, MBTET, "MESH_NODE_POSITIONS", 2);
376  CHKERR m_field.set_field_order(0, MBTRI, "MESH_NODE_POSITIONS", 2);
377  CHKERR m_field.set_field_order(0, MBEDGE, "MESH_NODE_POSITIONS", 2);
378  CHKERR m_field.set_field_order(0, MBVERTEX, "MESH_NODE_POSITIONS", 1);
379 
380  bool check_if_spatial_field_exist = m_field.check_field("DISPLACEMENT");
381  CHKERR m_field.add_field("DISPLACEMENT", H1, base, 3, MB_TAG_SPARSE,
382  MF_ZERO);
383  // add entities (by tets) to the field
384  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DISPLACEMENT");
385 
386  // set app. order
387  CHKERR m_field.set_field_order(0, MBTET, "DISPLACEMENT", disp_order);
388  CHKERR m_field.set_field_order(0, MBTRI, "DISPLACEMENT", disp_order);
389  CHKERR m_field.set_field_order(0, MBEDGE, "DISPLACEMENT", disp_order);
391  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", disp_order);
392  else
393  CHKERR m_field.set_field_order(0, MBVERTEX, "DISPLACEMENT", 1);
394 
395  // Add nodal force element
396  CHKERR MetaNeumannForces::addNeumannBCElements(m_field, "DISPLACEMENT");
397  CHKERR MetaEdgeForces::addElement(m_field, "DISPLACEMENT");
398  CHKERR MetaNodalForces::addElement(m_field, "DISPLACEMENT");
399  // Add fluid pressure finite elements
400  FluidPressure fluid_pressure_fe(m_field);
401  fluid_pressure_fe.addNeumannFluidPressureBCElements("DISPLACEMENT");
402  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS", fluid_pressure_fe.getLoopFe(),
403  false, false);
405  "DISPLACEMENT", PETSC_NULL, false, true);
406 
407  // Velocity
408  CHKERR m_field.add_field("VELOCITY", H1, base, 3, MB_TAG_SPARSE, MF_ZERO);
409  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "VELOCITY");
410 
411  CHKERR m_field.set_field_order(0, MBTET, "VELOCITY", vel_order);
412  CHKERR m_field.set_field_order(0, MBTRI, "VELOCITY", vel_order);
413  CHKERR m_field.set_field_order(0, MBEDGE, "VELOCITY", vel_order);
415  CHKERR m_field.set_field_order(0, MBVERTEX, "VELOCITY", vel_order);
416  else
417  CHKERR m_field.set_field_order(0, MBVERTEX, "VELOCITY", 1);
418 
419  CHKERR m_field.add_field("DOT_DISPLACEMENT", H1, base, 3, MB_TAG_SPARSE,
420  MF_ZERO);
421  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DOT_DISPLACEMENT");
422  CHKERR m_field.set_field_order(0, MBTET, "DOT_DISPLACEMENT", disp_order);
423  CHKERR m_field.set_field_order(0, MBTRI, "DOT_DISPLACEMENT", disp_order);
424  CHKERR m_field.set_field_order(0, MBEDGE, "DOT_DISPLACEMENT", disp_order);
426  CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_DISPLACEMENT",
427  disp_order);
428  else
429  CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_DISPLACEMENT", 1);
430 
431  CHKERR m_field.add_field("DOT_VELOCITY", H1, base, 3, MB_TAG_SPARSE,
432  MF_ZERO);
433  CHKERR m_field.add_ents_to_field_by_type(0, MBTET, "DOT_VELOCITY");
434  CHKERR m_field.set_field_order(0, MBTET, "DOT_VELOCITY", vel_order);
435  CHKERR m_field.set_field_order(0, MBTRI, "DOT_VELOCITY", vel_order);
436  CHKERR m_field.set_field_order(0, MBEDGE, "DOT_VELOCITY", vel_order);
438  CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_VELOCITY", disp_order);
439  else
440  CHKERR m_field.set_field_order(0, MBVERTEX, "DOT_VELOCITY", 1);
441 
442  // Set material model and mass element
443  NonlinearElasticElement elastic(m_field, 2);
444  ElasticMaterials elastic_materials(m_field);
445  CHKERR elastic_materials.setBlocks(elastic.setOfBlocks);
446  // NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<adouble>
447  // st_venant_kirchhoff_material_adouble;
448  // NonlinearElasticElement::FunctionsToCalculatePiolaKirchhoffI<double>
449  // st_venant_kirchhoff_material_double; CHKERR
450  // elastic.setBlocks(&st_venant_kirchhoff_material_double,&st_venant_kirchhoff_material_adouble);
451  CHKERR elastic.addElement("ELASTIC", "DISPLACEMENT");
452  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeRhs(), true, false,
453  false, false);
454  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeLhs(), true, false,
455  false, false);
456  CHKERR addHOOpsVol("MESH_NODE_POSITIONS", elastic.getLoopFeEnergy(), true,
457  false, false, false);
458  CHKERR elastic.setOperators("DISPLACEMENT", "MESH_NODE_POSITIONS", false,
459  true);
460 
461  // set mass element
462  ConvectiveMassElement inertia(m_field, 1);
463  // CHKERR inertia.setBlocks();
464  CHKERR elastic_materials.setBlocks(inertia.setOfBlocks);
465  CHKERR inertia.addConvectiveMassElement("MASS_ELEMENT", "VELOCITY",
466  "DISPLACEMENT");
467  CHKERR inertia.addVelocityElement("VELOCITY_ELEMENT", "VELOCITY",
468  "DISPLACEMENT");
469 
470  // Add possibility to load accelerogram
471  {
472  string name = "-my_accelerogram";
473  char time_file_name[255];
474  PetscBool flg;
475  CHKERR PetscOptionsGetString(PETSC_NULL, PETSC_NULL, name.c_str(),
476  time_file_name, 255, &flg);
477  if (flg == PETSC_TRUE) {
478  inertia.methodsOp.push_back(new TimeAccelerogram(name));
479  }
480  }
481 
482  // damper element
483  KelvinVoigtDamper damper(m_field);
484  CHKERR elastic_materials.setBlocks(damper.blockMaterialDataMap);
485  {
486  KelvinVoigtDamper::CommonData &common_data = damper.commonData;
487  common_data.spatialPositionName = "DISPLACEMENT";
488  common_data.spatialPositionNameDot = "DOT_DISPLACEMENT";
489  CHKERR m_field.add_finite_element("DAMPER", MF_ZERO);
491  "DISPLACEMENT");
493  "DISPLACEMENT");
495  "DISPLACEMENT");
497  "DOT_DISPLACEMENT");
498  if (m_field.check_field("MESH_NODE_POSITIONS")) {
500  "DAMPER", "MESH_NODE_POSITIONS");
501  }
502  std::map<int, KelvinVoigtDamper::BlockMaterialData>::iterator bit =
503  damper.blockMaterialDataMap.begin();
504  for (; bit != damper.blockMaterialDataMap.end(); bit++) {
505  bit->second.lInear = linear;
506  int id = bit->first;
507  KelvinVoigtDamper::BlockMaterialData &material_data = bit->second;
508  damper.constitutiveEquationMap.insert(
510  material_data));
511  CHKERR m_field.add_ents_to_finite_element_by_type(bit->second.tEts,
512  MBTET, "DAMPER");
513  }
514  CHKERR damper.setOperators(3);
515  }
516 
517  MonitorPostProc post_proc(m_field, elastic.setOfBlocks,
518  elastic.getLoopFeEnergy(),
519  inertia.getLoopFeEnergy());
520 
521  // elastic and mass element calculated in Kuu shell matrix problem. To
522  // calculate Mass element, velocity field is needed.
523  CHKERR m_field.modify_finite_element_add_field_data("ELASTIC", "VELOCITY");
525  "DOT_DISPLACEMENT");
527  "DOT_VELOCITY");
528 
529  // build field
530  CHKERR m_field.build_fields();
531  // CHKERR m_field.list_dofs_by_field_name("DISPLACEMENT");
532 
533  // 10 node tets
534  if (!check_if_spatial_field_exist) {
535  Projection10NodeCoordsOnField ent_method_material(m_field,
536  "MESH_NODE_POSITIONS");
537  CHKERR m_field.loop_dofs("MESH_NODE_POSITIONS", ent_method_material);
538  }
539 
540  // build finite elements
541  CHKERR m_field.build_finite_elements();
542  // build adjacencies
543  CHKERR m_field.build_adjacencies(bit_level0);
544 
545  // define problems
546  {
547  CHKERR m_field.add_problem("Kuu", MF_ZERO);
548  CHKERR m_field.modify_problem_add_finite_element("Kuu", "ELASTIC");
549  CHKERR m_field.modify_problem_add_finite_element("Kuu", "PRESSURE_FE");
550  CHKERR m_field.modify_problem_add_finite_element("Kuu", "FORCE_FE");
552  "FLUID_PRESSURE_FE");
553  CHKERR m_field.modify_problem_ref_level_add_bit("Kuu", bit_level0);
554 
555  ProblemsManager *prb_mng_ptr;
556  CHKERR m_field.getInterface(prb_mng_ptr);
557  if (is_partitioned) {
558  CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("Kuu", true);
559  CHKERR prb_mng_ptr->partitionFiniteElements("Kuu", true, 0,
560  pcomm->size());
561  } else {
562  CHKERR prb_mng_ptr->buildProblem("Kuu", true);
563  CHKERR prb_mng_ptr->partitionProblem("Kuu");
564  CHKERR prb_mng_ptr->partitionFiniteElements("Kuu");
565  }
566  CHKERR prb_mng_ptr->partitionGhostDofs("Kuu");
567  }
568 
569  CHKERR m_field.add_problem("DYNAMICS", MF_ZERO);
570  // set finite elements for problems
571  CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "ELASTIC");
572  CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "DAMPER");
573  CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "PRESSURE_FE");
574  CHKERR m_field.modify_problem_add_finite_element("DYNAMICS", "FORCE_FE");
575  CHKERR m_field.modify_problem_add_finite_element("DYNAMICS",
576  "FLUID_PRESSURE_FE");
577  CHKERR m_field.modify_problem_add_finite_element("DYNAMICS",
578  "MASS_ELEMENT");
579  CHKERR m_field.modify_problem_add_finite_element("DYNAMICS",
580  "VELOCITY_ELEMENT");
581  // set refinement level for problem
582  CHKERR m_field.modify_problem_ref_level_add_bit("DYNAMICS", bit_level0);
583 
584  ProblemsManager *prb_mng_ptr;
585  CHKERR m_field.getInterface(prb_mng_ptr);
586  if (is_partitioned) {
587  CHKERR prb_mng_ptr->buildProblemOnDistributedMesh("DYNAMICS", true);
588  CHKERR prb_mng_ptr->partitionFiniteElements("DYNAMICS", true, 0,
589  pcomm->size());
590  } else {
591  CHKERR prb_mng_ptr->buildProblem("DYNAMICS", true);
592  CHKERR prb_mng_ptr->partitionProblem("DYNAMICS");
593  CHKERR prb_mng_ptr->partitionFiniteElements("DYNAMICS");
594  }
595  CHKERR prb_mng_ptr->partitionGhostDofs("DYNAMICS");
596 
597  Vec F;
598  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("DYNAMICS", COL,
599  &F);
600  Vec D;
601  CHKERR VecDuplicate(F, &D);
602 
603  // create tS
604  TS ts;
605  CHKERR TSCreate(PETSC_COMM_WORLD, &ts);
606  CHKERR TSSetType(ts, TSBEULER);
607 
608  // shell matrix
609  ConvectiveMassElement::MatShellCtx *shellAij_ctx =
612  ->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("Kuu",
613  &shellAij_ctx->K);
614  CHKERR MatDuplicate(shellAij_ctx->K, MAT_DO_NOT_COPY_VALUES,
615  &shellAij_ctx->M);
616  CHKERR shellAij_ctx->iNit();
617  CHKERR m_field.getInterface<VecManager>()->vecScatterCreate(
618  D, "DYNAMICS", COL, shellAij_ctx->u, "Kuu", COL,
619  &shellAij_ctx->scatterU);
620  CHKERR m_field.getInterface<VecManager>()->vecScatterCreate(
621  D, "DYNAMICS", "VELOCITY", COL, shellAij_ctx->v, "Kuu", "DISPLACEMENT",
622  COL, &shellAij_ctx->scatterV);
623  Mat shell_Aij;
624  const Problem *problem_ptr;
625  CHKERR m_field.get_problem("DYNAMICS", &problem_ptr);
626  CHKERR MatCreateShell(
627  PETSC_COMM_WORLD, problem_ptr->getNbLocalDofsRow(),
628  problem_ptr->getNbLocalDofsCol(), problem_ptr->getNbDofsRow(),
629  problem_ptr->getNbDofsRow(), (void *)shellAij_ctx, &shell_Aij);
630  CHKERR MatShellSetOperation(shell_Aij, MATOP_MULT,
631  (void (*)(void))ConvectiveMassElement::MultOpA);
632  CHKERR MatShellSetOperation(
633  shell_Aij, MATOP_ZERO_ENTRIES,
634  (void (*)(void))ConvectiveMassElement::ZeroEntriesOp);
635  // blocked problem
636  ConvectiveMassElement::ShellMatrixElement shell_matrix_element(m_field);
637  DirichletDisplacementBc shell_dirichlet_bc(
638  m_field, "DISPLACEMENT", shellAij_ctx->barK, PETSC_NULL, PETSC_NULL);
639  DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT", PETSC_NULL,
640  D, F);
641  shell_matrix_element.problemName = "Kuu";
642  shell_matrix_element.shellMatCtx = shellAij_ctx;
643  shell_matrix_element.DirichletBcPtr = &shell_dirichlet_bc;
644  shell_matrix_element.loopK.push_back(
645  ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
646  "ELASTIC", &elastic.getLoopFeLhs()));
647  // damper
648  shell_matrix_element.loopK.push_back(
649  ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
650  "ELASTIC", &damper.feLhs));
651 
652  CHKERR inertia.addHOOpsVol();
653  CHKERR inertia.setShellMatrixMassOperators("VELOCITY", "DISPLACEMENT",
654  "MESH_NODE_POSITIONS", linear);
655  // element name "ELASTIC" is used, therefore M matrix is assembled as K
656  // matrix. This is added to M is shell matrix. M matrix is a derivative of
657  // inertia forces over spatial velocities
658  shell_matrix_element.loopM.push_back(
659  ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
660  "ELASTIC", &inertia.getLoopFeMassLhs()));
661  // this calculate derivatives of inertia forces over spatial positions and
662  // add this to shell K matrix
663  shell_matrix_element.loopAuxM.push_back(
664  ConvectiveMassElement::ShellMatrixElement::PairNameFEMethodPtr(
665  "ELASTIC", &inertia.getLoopFeMassAuxLhs()));
666 
667  // Element to calculate shell matrix residual
668  ConvectiveMassElement::ShellResidualElement shell_matrix_residual(m_field);
669  shell_matrix_residual.shellMatCtx = shellAij_ctx;
670 
671  // surface pressure
672  boost::ptr_map<std::string, NeumannForcesSurface> surface_forces;
673  {
674  string fe_name_str = "FORCE_FE";
675  surface_forces.insert(fe_name_str, new NeumannForcesSurface(m_field));
676  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS",
677  surface_forces.at(fe_name_str).getLoopFe(), false,
678  false);
680  NODESET | FORCESET, it)) {
681  CHKERR surface_forces.at(fe_name_str)
682  .addForce("DISPLACEMENT", PETSC_NULL, it->getMeshsetId(), true);
683  surface_forces.at(fe_name_str)
684  .methodsOp.push_back(new TimeForceScale());
685  }
686  }
687 
688  boost::ptr_map<std::string, NeumannForcesSurface> surface_pressure;
689  {
690  string fe_name_str = "PRESSURE_FE";
691  surface_pressure.insert(fe_name_str, new NeumannForcesSurface(m_field));
692  CHKERR addHOOpsFace3D("MESH_NODE_POSITIONS",
693  surface_pressure.at(fe_name_str).getLoopFe(), false,
694  false);
696  m_field, SIDESET | PRESSURESET, it)) {
697  CHKERR surface_pressure.at(fe_name_str)
698  .addPressure("DISPLACEMENT", PETSC_NULL, it->getMeshsetId(), true);
699  surface_pressure.at(fe_name_str)
700  .methodsOp.push_back(new TimeForceScale());
701  }
702  }
703 
704  // edge forces
705  boost::ptr_map<std::string, EdgeForce> edge_forces;
706  {
707  string fe_name_str = "FORCE_FE";
708  edge_forces.insert(fe_name_str, new EdgeForce(m_field));
710  NODESET | FORCESET, it)) {
711  CHKERR edge_forces.at(fe_name_str)
712  .addForce("DISPLACEMENT", PETSC_NULL, it->getMeshsetId(), true);
713  edge_forces.at(fe_name_str).methodsOp.push_back(new TimeForceScale());
714  }
715  }
716 
717  // nodal forces
718  boost::ptr_map<std::string, NodalForce> nodal_forces;
719  {
720  string fe_name_str = "FORCE_FE";
721  nodal_forces.insert(fe_name_str, new NodalForce(m_field));
723  NODESET | FORCESET, it)) {
724  CHKERR nodal_forces.at(fe_name_str)
725  .addForce("DISPLACEMENT", F, it->getMeshsetId(), true);
726  nodal_forces.at(fe_name_str).methodsOp.push_back(new TimeForceScale());
727  }
728  }
729 
730  MonitorRestart monitor_restart(m_field, ts);
731  ConvectiveMassElement::UpdateAndControl update_and_control(
732  m_field, ts, "VELOCITY", "DISPLACEMENT");
733 
734  // TS
735  TsCtx ts_ctx(m_field, "DYNAMICS");
736 
737  // right hand side
738  // preprocess
739  ts_ctx.getPreProcessIFunction().push_back(&update_and_control);
740  ts_ctx.getPreProcessIFunction().push_back(&my_dirichlet_bc);
741 
742  // fe looops
743  auto &loops_to_do_Rhs = ts_ctx.getLoopsIFunction();
744 
745  auto add_static_rhs = [&](auto &loops_to_do_Rhs) {
747  loops_to_do_Rhs.push_back(
748  PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeRhs()));
749  for (auto fit = surface_forces.begin(); fit != surface_forces.end();
750  fit++) {
751  loops_to_do_Rhs.push_back(
752  PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
753  }
754  for (auto fit = surface_pressure.begin(); fit != surface_pressure.end();
755  fit++) {
756  loops_to_do_Rhs.push_back(
757  PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
758  }
759  for (auto fit = edge_forces.begin(); fit != edge_forces.end(); fit++) {
760  loops_to_do_Rhs.push_back(
761  PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
762  }
763  for (auto fit = nodal_forces.begin(); fit != nodal_forces.end(); fit++) {
764  loops_to_do_Rhs.push_back(
765  PairNameFEMethodPtr(fit->first, &fit->second->getLoopFe()));
766  }
767  loops_to_do_Rhs.push_back(PairNameFEMethodPtr(
768  "FLUID_PRESSURE_FE", &fluid_pressure_fe.getLoopFe()));
770  };
771 
772  CHKERR add_static_rhs(loops_to_do_Rhs);
773 
774  loops_to_do_Rhs.push_back(PairNameFEMethodPtr("DAMPER", &damper.feRhs));
775  loops_to_do_Rhs.push_back(
776  PairNameFEMethodPtr("MASS_ELEMENT", &inertia.getLoopFeMassRhs()));
777 
778  // preporcess
779  // calculate residual for velocities
780  ts_ctx.getPreProcessIFunction().push_back(&shell_matrix_residual);
781  // postprocess
782  ts_ctx.getPostProcessIFunction().push_back(&my_dirichlet_bc);
783 
784  // left hand side
785  // preprocess
786  ts_ctx.getPreProcessIJacobian().push_back(&update_and_control);
787  ts_ctx.getPreProcessIJacobian().push_back(&shell_matrix_element);
788  ts_ctx.getPostProcessIJacobian().push_back(&update_and_control);
789  // monitor
790  TsCtx::FEMethodsSequence &loopsMonitor =
791  ts_ctx.getLoopsMonitor();
792  loopsMonitor.push_back(
793  TsCtx::PairNameFEMethodPtr("MASS_ELEMENT", &post_proc));
794  loopsMonitor.push_back(
795  TsCtx::PairNameFEMethodPtr("MASS_ELEMENT", &monitor_restart));
796 
797  CHKERR TSSetIFunction(ts, F, TsSetIFunction, &ts_ctx);
798  CHKERR TSSetIJacobian(ts, shell_Aij, shell_Aij, TsSetIJacobian, &ts_ctx);
799 
800  CHKERR TSMonitorSet(ts, TsMonitorSet, &ts_ctx, PETSC_NULL);
801 
802  double ftime = 1;
803  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
804  CHKERR TSSetSolution(ts, D);
805  CHKERR TSSetFromOptions(ts);
806  // shell matrix pre-conditioner
807  SNES snes;
808  CHKERR TSGetSNES(ts, &snes);
809  // CHKERR SNESSetFromOptions(snes);
810  KSP ksp;
811  CHKERR SNESGetKSP(snes, &ksp);
812  CHKERR KSPSetFromOptions(ksp);
813  PC pc;
814  CHKERR KSPGetPC(ksp, &pc);
815  CHKERR PCSetType(pc, PCSHELL);
816  ConvectiveMassElement::PCShellCtx pc_shell_ctx(shell_Aij);
817  CHKERR PCShellSetContext(pc, (void *)&pc_shell_ctx);
818  CHKERR PCShellSetApply(pc, ConvectiveMassElement::PCShellApplyOp);
819  CHKERR PCShellSetSetUp(pc, ConvectiveMassElement::PCShellSetUpOp);
820  CHKERR PCShellSetDestroy(pc, ConvectiveMassElement::PCShellDestroy);
821 
822  CHKERR VecZeroEntries(D);
823  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
824  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
825  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
826  "DYNAMICS", COL, D, INSERT_VALUES, SCATTER_REVERSE);
827 
828  // Solve problem at time Zero
829  if (is_solve_at_time_zero) {
830 
831  Mat Aij = shellAij_ctx->K;
832  Vec F;
833  CHKERR m_field.getInterface<VecManager>()->vecCreateGhost("Kuu", COL, &F);
834  Vec D;
835  CHKERR VecDuplicate(F, &D);
836 
837  // Set vector for Kuu problem from the mesh data
838  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
839  "Kuu", COL, D, INSERT_VALUES, SCATTER_FORWARD);
840  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
841  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
842 
843  SnesCtx snes_ctx(m_field, "Kuu");
844 
845  SNES snes;
846  CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
847  CHKERR SNESSetApplicationContext(snes, &snes_ctx);
848  CHKERR SNESSetFunction(snes, F, SnesRhs, &snes_ctx);
849  CHKERR SNESSetJacobian(snes, Aij, Aij, SnesMat, &snes_ctx);
850  CHKERR SNESSetFromOptions(snes);
851 
852  DirichletDisplacementBc my_dirichlet_bc(m_field, "DISPLACEMENT",
853  PETSC_NULL, D, F);
854 
855  SnesCtx::FEMethodsSequence &loops_to_do_Rhs =
856  snes_ctx.get_loops_to_do_Rhs();
857  snes_ctx.get_preProcess_to_do_Rhs().push_back(&my_dirichlet_bc);
858  fluid_pressure_fe.getLoopFe().ts_t = 0;
859  CHKERR add_static_rhs(loops_to_do_Rhs);
860  snes_ctx.get_postProcess_to_do_Rhs().push_back(&my_dirichlet_bc);
861 
862  SnesCtx::FEMethodsSequence &loops_to_do_Mat =
863  snes_ctx.get_loops_to_do_Mat();
864  snes_ctx.get_preProcess_to_do_Mat().push_back(&my_dirichlet_bc);
865  loops_to_do_Mat.push_back(
866  SnesCtx::PairNameFEMethodPtr("ELASTIC", &elastic.getLoopFeLhs()));
867  snes_ctx.get_postProcess_to_do_Mat().push_back(&my_dirichlet_bc);
868 
869  CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0, "VELOCITY");
870  CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0,
871  "DOT_DISPLACEMENT");
872  CHKERR m_field.getInterface<FieldBlas>()->fieldScale(0, "DOT_VELOCITY");
873 
874  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
875  "Kuu", COL, D, INSERT_VALUES, SCATTER_FORWARD);
876 
877  CHKERR SNESSolve(snes, PETSC_NULL, D);
878  int its;
879  CHKERR SNESGetIterationNumber(snes, &its);
880  MOFEM_LOG_C("DYNAMIC", Sev::inform, "number of Newton iterations = %d\n",
881  its);
882 
883  // Set data on the mesh
884  CHKERR m_field.getInterface<VecManager>()->setGlobalGhostVector(
885  "Kuu", COL, D, INSERT_VALUES, SCATTER_REVERSE);
886 
887  CHKERR VecDestroy(&F);
888  CHKERR VecDestroy(&D);
889  CHKERR SNESDestroy(&snes);
890  }
891 
892  if (is_solve_at_time_zero) {
893  CHKERR m_field.getInterface<VecManager>()->setLocalGhostVector(
894  "DYNAMICS", COL, D, INSERT_VALUES, SCATTER_FORWARD);
895  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
896  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
897  CHKERR TSSetSolution(ts, D);
898  }
899 
900 #if PETSC_VERSION_GE(3, 7, 0)
901  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
902 #endif
903  CHKERR TSSolve(ts, D);
904  CHKERR TSGetTime(ts, &ftime);
905 
906  PetscInt steps, snesfails, rejects, nonlinits, linits;
907  CHKERR TSGetTimeStepNumber(ts, &steps);
908  CHKERR TSGetSNESFailures(ts, &snesfails);
909  CHKERR TSGetStepRejections(ts, &rejects);
910  CHKERR TSGetSNESIterations(ts, &nonlinits);
911  CHKERR TSGetKSPIterations(ts, &linits);
912  MOFEM_LOG_C("DYNAMIC", Sev::inform,
913  "steps %d (%d rejected, %D SNES fails), ftime %g, nonlinits "
914  "%d, linits %D\n",
915  steps, rejects, snesfails, ftime, nonlinits, linits);
916  CHKERR TSDestroy(&ts);
917 
918  CHKERR VecDestroy(&F);
919  CHKERR VecDestroy(&D);
920  CHKERR MatDestroy(&shellAij_ctx->K);
921  CHKERR MatDestroy(&shellAij_ctx->M);
922  CHKERR VecScatterDestroy(&shellAij_ctx->scatterU);
923  CHKERR VecScatterDestroy(&shellAij_ctx->scatterV);
924  CHKERR MatDestroy(&shell_Aij);
925  delete shellAij_ctx;
926  }
927  CATCH_ERRORS;
928 
930 
931  return 0;
932 }
const std::string default_options
std::string param_file
Elastic materials.
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:314
@ COL
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ MF_ZERO
Definition: definitions.h:111
FieldApproximationBase
approximation base
Definition: definitions.h:71
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ AINSWORTH_LOBATTO_BASE
Definition: definitions.h:75
@ NOBASE
Definition: definitions.h:72
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:608
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#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
@ PRESSURESET
Definition: definitions.h:165
@ FORCESET
Definition: definitions.h:164
@ NODESET
Definition: definitions.h:159
@ SIDESET
Definition: definitions.h:160
#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
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
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
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode add_ents_to_finite_element_by_type(const EntityHandle entities, const EntityType type, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
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.
virtual MoFEMErrorCode add_ents_to_field_by_type(const Range &ents, const EntityType type, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual bool check_field(const std::string &name) const =0
check if field is in database
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
MoFEMErrorCode addFieldValuesGradientPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2 or H1 field gradient.
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
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.
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.
#define _IT_CUBITMESHSETS_BY_BCDATA_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet in a moFEM field.
MoFEMErrorCode buildProblemOnDistributedMesh(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures, assuming that mesh is distributed (collective)
MoFEMErrorCode partitionGhostDofs(const std::string name, int verb=VERBOSE)
determine ghost nodes
MoFEMErrorCode buildProblem(const std::string name, const bool square_matrix, int verb=VERBOSE)
build problem data structures
MoFEMErrorCode partitionProblem(const std::string name, int verb=VERBOSE)
partition problem dofs (collective)
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
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
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
auto bit
set bit
char mesh_file_name[255]
const double T
const FTensor::Tensor2< T, Dim, Dim > Vec
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
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
PetscErrorCode TsSetIJacobian(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Set function evaluating jacobina in TS solver.
Definition: TsCtx.cpp:140
PetscErrorCode TsMonitorSet(TS ts, PetscInt step, PetscReal t, Vec u, void *ctx)
Set monitor for TS solver.
Definition: TsCtx.cpp:229
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:148
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
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:39
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode addHOOpsFace3D(const std::string field, E &e, bool hcurl, bool hdiv)
PetscErrorCode TsSetIFunction(TS ts, PetscReal t, Vec u, Vec u_t, Vec F, void *ctx)
Set IFunction for TS solver.
Definition: TsCtx.cpp:37
CoreTmp< 0 > Core
Definition: Core.hpp:1096
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
MoFEMErrorCode addHOOpsVol(const std::string field, E &e, bool h1, bool hcurl, bool hdiv, bool l2)
const double D
diffusivity
int main(int argc, char *argv[])
static char help[]
MatShellCtx * shellMatCtx
pointer to shell matrix
structure grouping operators and data used for calculation of mass (convective) element \ nonlinear_e...
static MoFEMErrorCode MultOpA(Mat A, Vec x, Vec f)
Mult operator for shell matrix.
static MoFEMErrorCode ZeroEntriesOp(Mat A)
MyVolumeFE & getLoopFeMassLhs()
get lhs volume element
MoFEMErrorCode setShellMatrixMassOperators(string velocity_field_name, string spatial_position_field_name, string material_position_field_name="MESH_NODE_POSITIONS", bool linear=false)
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())
static MoFEMErrorCode PCShellDestroy(PC pc)
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())
static MoFEMErrorCode PCShellApplyOp(PC pc, Vec f, Vec x)
apply pre-conditioner for shell matrix
MyVolumeFE & getLoopFeMassAuxLhs()
get lhs volume element for Kuu shell matrix
boost::ptr_vector< MethodForForceScaling > methodsOp
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
static MoFEMErrorCode PCShellSetUpOp(PC pc)
MyVolumeFE & getLoopFeMassRhs()
get rhs volume element
MyVolumeFE & getLoopFeEnergy()
get kinetic energy element
Set Dirichlet boundary conditions on displacements.
Definition: DirichletBC.hpp:69
Force on edges and lines.
Definition: EdgeForce.hpp:25
Manage setting parameters and constitutive equations for nonlinear/linear elastic materials.
Fluid pressure forces.
MyTriangleFE & getLoopFe()
MoFEMErrorCode setNeumannFluidPressureFiniteElementOperators(string field_name, Vec F, bool allow_negative_pressure=true, bool ho_geometry=false)
MoFEMErrorCode addNeumannFluidPressureBCElements(const std::string field_name, const std::string mesh_nodals_positions="MESH_NODE_POSITIONS")
Common data for nonlinear_elastic_elem model.
Implementation of Kelvin Voigt Damper.
ConstitutiveEquationMap constitutiveEquationMap
std::map< int, BlockMaterialData > blockMaterialDataMap
MoFEMErrorCode setOperators(const int tag)
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 addElement(MoFEM::Interface &m_field, const std::string field_name, Range *intersect_ptr=NULL)
Add element taking information from NODESET.
Definition: NodalForce.hpp:104
Managing BitRefLevels.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
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
Add field.
virtual moab::Interface & get_moab()=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
Basic algebra on fields.
Definition: FieldBlas.hpp:31
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
Matrix manager is used to build and partition problems.
keeps basic data about problem
DofIdx getNbLocalDofsRow() const
DofIdx getNbDofsRow() const
DofIdx getNbLocalDofsCol() const
Problem manager is used to build and partition problems.
Projection of edge entities with one mid-node on hierarchical basis.
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:27
BasicMethodsSequence & get_postProcess_to_do_Rhs()
Definition: SnesCtx.hpp:95
BasicMethodsSequence & get_preProcess_to_do_Mat()
Definition: SnesCtx.hpp:100
BasicMethodsSequence & get_preProcess_to_do_Rhs()
Definition: SnesCtx.hpp:86
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: SnesCtx.hpp:41
FEMethodsSequence & get_loops_to_do_Rhs()
Definition: SnesCtx.hpp:77
FEMethodsSequence & get_loops_to_do_Mat()
Definition: SnesCtx.hpp:72
BasicMethodsSequence & get_postProcess_to_do_Mat()
Definition: SnesCtx.hpp:109
Interface for Time Stepping (TS) solver.
Definition: TsCtx.hpp:27
FEMethodsSequence & getLoopsIFunction()
Get the loops to do IFunction object.
Definition: TsCtx.hpp:81
MoFEM::FEMethodsSequence FEMethodsSequence
Definition: TsCtx.hpp:36
FEMethodsSequence & getLoopsMonitor()
Get the loops to do Monitor object.
Definition: TsCtx.hpp:120
BasicMethodsSequence & getPostProcessIFunction()
Get the postProcess to do IFunction object.
Definition: TsCtx.hpp:134
BasicMethodsSequence & getPreProcessIFunction()
Get the preProcess to do IFunction object.
Definition: TsCtx.hpp:127
BasicMethodsSequence & getPostProcessIJacobian()
Get the postProcess to do IJacobian object.
Definition: TsCtx.hpp:150
BasicMethodsSequence & getPreProcessIJacobian()
Get the preProcess to do IJacobian object.
Definition: TsCtx.hpp:143
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Vector manager is used to create vectors \mofem_vectors.
Definition: VecManager.hpp:33
NonlinearElasticElement::MyVolumeFE & feElasticEnergy
calculate elastic energy
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEM::Interface & mField
std::map< int, NonlinearElasticElement::BlockData > & setOfBlocks
ConvectiveMassElement::MyVolumeFE & feKineticEnergy
calculate elastic energy
bool iNit
int pRT
MonitorPostProc(MoFEM::Interface &m_field, std::map< int, NonlinearElasticElement::BlockData > &set_of_blocks, NonlinearElasticElement::MyVolumeFE &fe_elastic_energy, ConvectiveMassElement::MyVolumeFE &fe_kinetic_energy)
PostProcVolumeOnRefinedMesh postProc
int * step
MonitorRestart(MoFEM::Interface &m_field, TS ts)
MoFEMErrorCode preProcess()
function is run at the beginning of loop
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
MoFEM::Interface & mField
Finite element and operators to apply force/pressures applied to surfaces.
Force applied to nodes.
Definition: NodalForce.hpp:25
structure grouping operators and data used for calculation of nonlinear elastic element
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)
std::map< int, BlockData > setOfBlocks
maps block set id with appropriate BlockData
MyVolumeFE & getLoopFeEnergy()
get energy fe
MyVolumeFE & getLoopFeLhs()
get lhs volume element
MyVolumeFE & getLoopFeRhs()
get rhs volume element
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.
std::vector< EntityHandle > mapGaussPts
moab::Interface & postProcMesh
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
CommonData commonData
Post processing.
Force scale operator for reading two columns.