v0.9.1
std_reac_diff.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
3 #include <rd_stdOperators.hpp>
4 
5 using namespace MoFEM;
6 using namespace StdRDOperators;
7 
8 static char help[] = "...\n\n";
9 
10 // struct RhsU {
11 // double operator()(const double u, const double v) const {
12 // return c * u * (u - alpha) * (1.0 - u) - u * v;
13 // }
14 // };
15 
16 // struct RhsV {
17 // double operator()(const double u, const double v) const {
18 // return (gma + mu1 * v / (mu2 + u)) * (-v - c * u * (u - b - 1.0));
19 // }
20 // };
21 
22 struct RDProblem {
23 public:
24  RDProblem(moab::Core &mb_instance, MoFEM::Core &core, const int order,
25  const int n_species)
26  : moab(mb_instance)
27  , m_field(core)
28  , order(order)
29  , nb_species(n_species)
30  , cOmm(m_field.get_comm())
31  , rAnk(m_field.get_comm_rank())
32  {
33  vol_ele_slow_rhs = boost::shared_ptr<Ele>(new Ele(m_field));
34  vol_ele_stiff_rhs = boost::shared_ptr<Ele>(new Ele(m_field));
35  vol_ele_stiff_lhs = boost::shared_ptr<Ele>(new Ele(m_field));
36 
37  boundary_ele_rhs = boost::shared_ptr<BoundaryEle>(new BoundaryEle(m_field));
38 
39  vol_mass_ele = boost::shared_ptr<Ele>(new Ele(m_field));
40 
41  post_proc = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
42  new PostProcFaceOnRefinedMesh(m_field));
43 
44  data.resize(nb_species);
45  values_ptr.resize(nb_species);
46  grads_ptr.resize(nb_species);
47  dots_ptr.resize(nb_species);
48  inner_surface.resize(nb_species);
49 
50  for (int i = 0; i < nb_species; ++i) {
51  data[i] = boost::shared_ptr<PreviousData>(new PreviousData());
52  grads_ptr[i] = boost::shared_ptr<MatrixDouble>(data[i], &data[i]->grads);
53  values_ptr[i] =
54  boost::shared_ptr<VectorDouble>(data[i], &data[i]->values);
55  dots_ptr[i] =
56  boost::shared_ptr<VectorDouble>(data[i], &data[i]->dot_values);
57  }
58  }
59 
60  // RDProblem(const int order) : order(order){}
61  MoFEMErrorCode run_analysis();
62 
63 private:
64  MoFEMErrorCode setup_system();
65  MoFEMErrorCode add_fe(std::string field_name);
66  MoFEMErrorCode set_blockData(std::map<int, BlockData> &block_data_map);
67 
68  MoFEMErrorCode set_initial_values(std::string field_name, int block_id,
69  Range &surface, double &init_val);
70 
71  MoFEMErrorCode update_slow_rhs(std::string mass_fiedl,
72  boost::shared_ptr<VectorDouble> &mass_ptr);
73 
74  MoFEMErrorCode push_slow_rhs(std::string mass_name,
75  boost::shared_ptr<PreviousData> &data);
76 
77  MoFEMErrorCode push_mass_ele(std::string field_name);
78 
79  MoFEMErrorCode resolve_slow_rhs();
80 
81  MoFEMErrorCode update_vol_fe(boost::shared_ptr<Ele> &vol_ele,
82  boost::shared_ptr<PreviousData> &data);
83  MoFEMErrorCode update_stiff_rhs(std::string field_name,
84  boost::shared_ptr<VectorDouble> &values_ptr,
85  boost::shared_ptr<MatrixDouble> &grads_ptr,
86  boost::shared_ptr<VectorDouble> &dots_ptr);
87 
88  MoFEMErrorCode push_stiff_rhs(std::string field_name,
89  boost::shared_ptr<PreviousData> &data,
90  std::map<int, BlockData> &block_map);
91 
92  MoFEMErrorCode update_stiff_lhs(std::string field_name,
93  boost::shared_ptr<VectorDouble> &values_ptr,
94  boost::shared_ptr<MatrixDouble> &grads_ptr);
95  MoFEMErrorCode push_stiff_lhs(std::string field_name,
96  boost::shared_ptr<PreviousData> &data,
97  std::map<int, BlockData> &block_map);
98 
99  MoFEMErrorCode set_integration_rule();
100 
101  MoFEMErrorCode set_fe_in_loop();
102  MoFEMErrorCode post_proc_fields(std::string field_name);
103  MoFEMErrorCode output_result();
104  MoFEMErrorCode solve();
105 
106  MoFEM::Interface &m_field;
107  Simple *simple_interface;
108 
110 
115 
116  Range bdry_ents;
117 
118  std::vector<Range> inner_surface;
120 
121  double global_error;
122 
123  MPI_Comm cOmm;
124  const int rAnk;
125 
126  int order;
127  int nb_species;
128 
129  std::map<int, BlockData> material_blocks;
130 
131  boost::shared_ptr<Ele> vol_ele_slow_rhs;
132  boost::shared_ptr<Ele> vol_ele_stiff_rhs;
133  boost::shared_ptr<Ele> vol_ele_stiff_lhs;
134  boost::shared_ptr<BoundaryEle> boundary_ele_rhs;
135 
136  boost::shared_ptr<Ele> vol_mass_ele;
137 
138  boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc;
139  boost::shared_ptr<Monitor> monitor_ptr;
140 
141  std::vector<boost::shared_ptr<PreviousData>> data;
142 
143  std::vector<boost::shared_ptr<MatrixDouble>> grads_ptr;
144  std::vector<boost::shared_ptr<VectorDouble>> values_ptr;
145  std::vector<boost::shared_ptr<VectorDouble>> dots_ptr;
146 
147  boost::shared_ptr<ForcesAndSourcesCore> null;
148 };
149 
150 const double ramp_t = 1.0;
151 const double sml = 0.0;
152 const double T = 2.0 * M_PI;
153 
154 struct ExactFunction {
155  double operator()(const double x, const double y, const double t) const {
156  double g = sin(T * x) * sin(T * y);
157  double val = 0;
158  if (x > -sml) {
159  val = 1.0 * g;
160  } else {
161  val = g;
162  }
163  if (t <= ramp_t) {
164  return val * t;
165  } else {
166  return val * ramp_t;
167  }
168  }
169 };
170 
171 struct ExactFunctionGrad {
172  FTensor::Tensor1<double, 3> operator()(const double x, const double y,
173  const double t) const {
175  double mx = -T * cos(T * x) * sin(T * y);
176  double my = -T * sin(T * x) * cos(T * y);
177  double hx, hy;
178  if (x > -sml) {
179  hx = 1.0 * mx;
180  hy = 1.0 * my;
181  } else {
182  hx = mx;
183  hy = my;
184  }
185  if (t <= ramp_t) {
186  grad(0) = hx * t;
187  grad(1) = hy * t;
188  } else {
189  grad(0) = hx * ramp_t;
190  grad(1) = hy * ramp_t;
191  }
192  grad(2) = 0.0;
193  return grad;
194  }
195 };
196 
197 struct ExactFunctionLap {
198  double operator()(const double x, const double y, const double t) const {
199  double glap = -2.0 * pow(T, 2) * sin(T * x) * sin(T * y);
200  double lap;
201  if (x > -sml) {
202  lap = 1.0 * glap;
203  } else {
204  lap = glap;
205  }
206  if (t <= ramp_t) {
207  return lap * t;
208  } else {
209  return lap * ramp_t;
210  }
211  }
212 };
213 
214 struct ExactFunctionDot {
215  double operator()(const double x, const double y, const double t) const {
216  double gdot = sin(T * x) * sin(T * y);
217  double dot;
218  if (x > -sml) {
219  dot = 1.0 * gdot;
220  } else {
221  dot = gdot;
222  }
223  if (t <= ramp_t) {
224  return dot;
225  } else {
226  return 0;
227  }
228  }
229 };
230 
231 // struct ReactionTerms {
232 // void operator()(const double u1, const double u2, const double u3,
233 // BlockData &coeffs, VectorDouble &vec){
234 // vec.resize(3, false);
235 // vec.clear();
236 
237 // vec[0] = coeffs.rate[0] * u1 * (1.0 - coeffs.coef(0, 0) * u1
238 // - coeffs.coef(0, 1) * u2
239 // - coeffs.coef(0, 2) * u3);
240 
241 // vec[1] = coeffs.rate[1] * u1 * (1.0 - coeffs.coef(1, 0) * u1
242 // - coeffs.coef(1, 1) * u2
243 // - coeffs.coef(1, 2) * u3);
244 // vec[2] = coeffs.rate[2] * (u1 * u2 + u3) * (1.0 - coeffs.coef(0, 0) * u1
245 // - coeffs.coef(0, 1) * u2
246 // - coeffs.coef(0, 2) * u3);
247 // }
248 
249 // };
252  CHKERR m_field.getInterface(simple_interface);
253  CHKERR simple_interface->getOptions();
254  CHKERR simple_interface->loadFile();
256 }
257 MoFEMErrorCode RDProblem::add_fe(std::string field_name) {
259  CHKERR simple_interface->addDomainField(field_name, H1,
261  CHKERR simple_interface->addDataField("ERROR", L2, AINSWORTH_LEGENDRE_BASE,
262  1);
263  // CHKERR simple_interface->addBoundaryField(field_name, H1,
264  // AINSWORTH_LEGENDRE_BASE, 1);
265  CHKERR simple_interface->setFieldOrder(field_name, order);
266  CHKERR simple_interface->setFieldOrder("ERROR", 0);
267 
269 }
270 
271 MoFEMErrorCode RDProblem::set_blockData(std::map<int, BlockData> &block_map) {
274  string name = it->getName();
275  const int id = it->getMeshsetId();
276  if (name.compare(0, 14, "REGION1") == 0) {
277  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
278  id, BLOCKSET, 2, block_map[id].block_ents, true);
279  block_map[id].B0 = 1e-3;
280  block_map[id].block_id = id;
281  } else if (name.compare(0, 14, "REGION2") == 0) {
282  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
283  id, BLOCKSET, 2, block_map[id].block_ents, true);
284  block_map[id].B0 = 1e-2;
285  block_map[id].block_id = id;
286  } else if (name.compare(0, 14, "REGION3") == 0) {
287  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
288  id, BLOCKSET, 2, block_map[id].block_ents, true);
289  block_map[id].B0 = 5e-2;
290  block_map[id].block_id = id;
291  } else if (name.compare(0, 14, "REGION4") == 0) {
292  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
293  id, BLOCKSET, 2, block_map[id].block_ents, true);
294  block_map[id].B0 = 1e-1;
295  block_map[id].block_id = id;
296  } else if (name.compare(0, 14, "REGION5") == 0) {
297  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
298  id, BLOCKSET, 2, block_map[id].block_ents, true);
299  block_map[id].B0 = 5e-1;
300  block_map[id].block_id = id;
301  }
302  }
304 }
305 
307  int block_id, Range &surface,
308  double &init_val) {
310  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(block_id,
311  BLOCKSET)) {
312  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
313  block_id, BLOCKSET, 2, surface, true);
314  }
315  if (!surface.empty()) {
316  Range surface_verts;
317  CHKERR moab.get_connectivity(surface, surface_verts, false);
318  CHKERR m_field.getInterface<FieldBlas>()->setField(
319  init_val, MBVERTEX, surface_verts, field_name);
320  }
321 
323 }
324 
326 RDProblem::update_slow_rhs(std::string mass_field,
327  boost::shared_ptr<VectorDouble> &values_ptr) {
329  vol_ele_slow_rhs->getOpPtrVector().push_back(
330  new OpCalculateScalarFieldValues(mass_field, values_ptr));
332 }
333 
335  boost::shared_ptr<PreviousData> &data) {
337 
338  vol_ele_slow_rhs->getOpPtrVector().push_back(
339  new OpAssembleSlowRhs(field_name, data, ExactFunction(),
342 }
343 
344 MoFEMErrorCode RDProblem::update_vol_fe(boost::shared_ptr<Ele> &vol_ele,
345  boost::shared_ptr<PreviousData> &data) {
347 
348  vol_ele->getOpPtrVector().push_back(
349  new OpCalculateInvJacForFace(data->invJac));
350  vol_ele->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(data->invJac));
351 
353 }
354 
356 RDProblem::update_stiff_rhs(std::string field_name,
357  boost::shared_ptr<VectorDouble> &values_ptr,
358  boost::shared_ptr<MatrixDouble> &grads_ptr,
359  boost::shared_ptr<VectorDouble> &dots_ptr) {
360 
362 
363  vol_ele_stiff_rhs->getOpPtrVector().push_back(
364  new OpCalculateScalarFieldValues(field_name, values_ptr));
365  vol_ele_stiff_rhs->getOpPtrVector().push_back(
366  new OpCalculateScalarValuesDot(field_name, dots_ptr));
367  vol_ele_stiff_rhs->getOpPtrVector().push_back(
368  new OpCalculateScalarFieldGradient<2>(field_name, grads_ptr));
370 }
372  boost::shared_ptr<PreviousData> &data,
373  std::map<int, BlockData> &block_map) {
375 
376  vol_ele_stiff_rhs->getOpPtrVector().push_back(
377  new OpAssembleStiffRhs<2>(field_name, data, block_map, ExactFunction(),
380 }
381 
383 RDProblem::update_stiff_lhs(std::string field_name,
384  boost::shared_ptr<VectorDouble> &values_ptr,
385  boost::shared_ptr<MatrixDouble> &grads_ptr) {
387  vol_ele_stiff_lhs->getOpPtrVector().push_back(
388  new OpCalculateScalarFieldValues(field_name, values_ptr));
389 
390  vol_ele_stiff_lhs->getOpPtrVector().push_back(
391  new OpCalculateScalarFieldGradient<2>(field_name, grads_ptr));
393 }
394 
396  boost::shared_ptr<PreviousData> &data,
397  std::map<int, BlockData> &block_map) {
399 
400  vol_ele_stiff_lhs->getOpPtrVector().push_back(
401  new OpAssembleStiffLhs<2>(field_name, data, block_map));
402 
404 }
405 
408  auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
409  vol_ele_slow_rhs->getRuleHook = vol_rule;
410 
411  vol_ele_stiff_rhs->getRuleHook = vol_rule;
412 
413  vol_ele_stiff_lhs->getRuleHook = vol_rule;
414  boundary_ele_rhs->getRuleHook = vol_rule;
416 }
417 
418 MoFEMErrorCode RDProblem::push_mass_ele(std::string field_name) {
420  vol_mass_ele->getOpPtrVector().push_back(
421  new OpAssembleMass(field_name, mass_matrix));
423 }
424 
427  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
428  vol_mass_ele);
429  CHKERR MatAssemblyBegin(mass_matrix, MAT_FINAL_ASSEMBLY);
430  CHKERR MatAssemblyEnd(mass_matrix, MAT_FINAL_ASSEMBLY);
431  // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
432  // M^-1G(t,u)
433  mass_Ksp = createKSP(m_field.get_comm());
434  CHKERR KSPSetOperators(mass_Ksp, mass_matrix, mass_matrix);
435  CHKERR KSPSetFromOptions(mass_Ksp);
436  CHKERR KSPSetUp(mass_Ksp);
438 }
439 
442  CHKERR TSSetType(ts, TSARKIMEX);
443  CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
444 
445  CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
446  vol_ele_stiff_lhs, null, null);
447 
448  CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
449  vol_ele_stiff_rhs, null, null);
450 
451  CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
452  vol_ele_slow_rhs, null, null);
454 }
455 
458  post_proc->addFieldValuesPostProc(field_name);
460 }
461 
464  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
465  monitor_ptr, null, null);
467 }
468 
471  // Create solution vector
474  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
475  // Solve problem
476  double ftime = 1;
477  CHKERR TSSetDM(ts, dm);
478  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
479  CHKERR TSSetSolution(ts, X);
480  CHKERR TSSetFromOptions(ts);
481  CHKERR TSSolve(ts, X);
483 }
484 
487  global_error = 0;
488  std::vector<std::string> mass_names(nb_species);
489 
490  for (int i = 0; i < nb_species; ++i) {
491  mass_names[i] = "MASS" + boost::lexical_cast<std::string>(i + 1);
492  }
493  CHKERR setup_system();
494 
495  for (int i = 0; i < nb_species; ++i) {
496  CHKERR add_fe(mass_names[i]);
497  }
498 
499  CHKERR simple_interface->setUp();
500 
501  // for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
502  // string name = it->getName();
503  // if (name.compare(0, 14, "ESSENTIAL") == 0) {
504  // CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
505  // bdry_ents, true);
506  // }
507  // }
508  Range surface;
509  CHKERR moab.get_entities_by_type(0, MBTRI, surface, false);
510  Skinner skin(&m_field.get_moab());
511  Range edges;
512  CHKERR skin.find_skin(0, surface, false, edges);
513  Range edges_part;
514  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
515  CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
516  PSTATUS_NOT, -1, &edges_part);
517  Range edges_verts;
518  CHKERR moab.get_connectivity(edges_part, edges_verts, false);
519 
520  CHKERR set_blockData(material_blocks);
521 
522  VectorDouble initVals;
523  initVals.resize(3, false);
524  initVals.clear();
525 
526  initVals[0] = 1.0;
527  initVals[1] = 3;
528  initVals[2] = 0.0;
529 
530  for (int i = 0; i < 1; ++i) {
531  CHKERR set_initial_values(mass_names[i], i + 2, inner_surface[i],
532  initVals[i]);
533  CHKERR update_slow_rhs(mass_names[i], values_ptr[i]);
534  }
535 
536  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(3, BLOCKSET)) {
537  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
538  3, BLOCKSET, 2, stimulation_surface, true);
539  }
540 
541  if (nb_species == 1) {
542  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
543  mass_names[0], data[0], data[0], data[0], material_blocks));
544  } else if (nb_species == 2) {
545  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
546  mass_names[0], data[0], data[1], data[1], material_blocks));
547  } else if (nb_species == 3) {
548  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
549  mass_names[0], data[0], data[1], data[2], material_blocks));
550  }
551 
552  for (int i = 0; i < nb_species; ++i) {
553  CHKERR push_slow_rhs(mass_names[i], data[i]);
554  // boundary_ele_rhs->getOpPtrVector().push_back(
555  // new OpAssembleNaturalBCRhs(mass_names[i], bdry_ents));
556  }
557 
558  dm = simple_interface->getDM();
559 
560  // cout << "essential : " << essential_bdry_ents.empty() << endl;
561 
562  auto solve_for_g = [&]() {
564  if (vol_ele_slow_rhs->vecAssembleSwitch) {
565  CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
566  SCATTER_REVERSE);
567  CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
568  SCATTER_REVERSE);
569  CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
570  CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
571  *vol_ele_slow_rhs->vecAssembleSwitch = false;
572  }
573  CHKERR KSPSolve(mass_Ksp, vol_ele_slow_rhs->ts_F, vol_ele_slow_rhs->ts_F);
575  };
576  // Add hook to the element to calculate g.
577  bdry_ents = unite(edges_verts, edges_part);
578 
579  // CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
580  // simple_interface->getProblemName(), mass_names[0], bdry_ents);
581 
582  CHKERR DMCreateMatrix_MoFEM(dm, mass_matrix);
583  CHKERR MatZeroEntries(mass_matrix);
584 
585  for (int i = 0; i < nb_species; ++i) {
586  CHKERR push_mass_ele(mass_names[i]);
587  }
588 
589  CHKERR resolve_slow_rhs();
590 
591  CHKERR update_vol_fe(vol_ele_stiff_rhs, data[0]);
592 
593  for (int i = 0; i < nb_species; ++i) {
594  CHKERR update_stiff_rhs(mass_names[i], values_ptr[i], grads_ptr[i],
595  dots_ptr[i]);
596  CHKERR push_stiff_rhs(mass_names[i], data[i], material_blocks);
597  }
598 
599  CHKERR update_vol_fe(vol_ele_stiff_lhs, data[0]);
600 
601  for (int i = 0; i < nb_species; ++i) {
602  CHKERR update_stiff_lhs(mass_names[i], values_ptr[i], grads_ptr[i]);
603  CHKERR push_stiff_lhs(mass_names[i], data[i],
604  material_blocks); // nb_species times
605  }
606  // vol_ele_stiff_lhs->getOpPtrVector().push_back(
607  // new OpError(ExactFunction(), ExactFunctionLap(), ExactFunctionGrad(),
608  // data[0], material_blocks, global_error));
609 
610  CHKERR set_integration_rule();
611 
612  ts = createTS(m_field.get_comm());
613 
614  vol_ele_slow_rhs->postProcessHook = solve_for_g;
615 
616  CHKERR set_fe_in_loop();
617 
618  post_proc->generateReferenceElementMesh(); // only once
619 
620  for (int i = 0; i < nb_species; ++i) {
621  CHKERR post_proc_fields(mass_names[i]);
622  }
623  post_proc->addFieldValuesPostProc("ERROR");
624 
625  monitor_ptr = boost::shared_ptr<Monitor>(
626  new Monitor(cOmm, rAnk, dm, post_proc, global_error)); // nb_species times
627  CHKERR output_result(); // only once
628 
629  CHKERR solve(); // only once
630 
632 }
633 
634 int main(int argc, char *argv[]) {
635  const char param_file[] = "param_file.petsc";
636  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
637  try {
638  moab::Core mb_instance;
639  moab::Interface &moab = mb_instance;
640  MoFEM::Core core(moab);
641  DMType dm_name = "DMMOFEM";
642  CHKERR DMRegister_MoFEM(dm_name);
643 
644  int order = 1;
645  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
646  int nb_species = 1;
647  RDProblem reac_diff_problem(mb_instance, core, order, nb_species);
648  CHKERR reac_diff_problem.run_analysis();
649  }
650  CATCH_ERRORS;
652  return 0;
653 }
FaceElementForcesAndSourcesCore Ele
MoFEMErrorCode resolve_slow_rhs()
std::vector< boost::shared_ptr< PreviousData > > data
RDProblem(moab::Core &mb_instance, MoFEM::Core &core, const int order, const int n_species)
boost::shared_ptr< Ele > vol_ele_stiff_lhs
Deprecated interface functions.
MoFEMErrorCode add_fe()
Range stimulation_surface
std::vector< boost::shared_ptr< MatrixDouble > > grads_ptr
int main(int argc, char *argv[])
MoFEMErrorCode push_mass_ele(std::string field_name)
MoFEMErrorCode update_stiff_lhs()
MoFEMErrorCode set_fe_in_loop()
MoFEMErrorCode setup_system()
auto createTS
Definition: AuxPETSc.hpp:275
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:445
double operator()(const double x, const double y, const double t) const
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
Interface for managing meshsets containing materials and boundary conditions.
boost::shared_ptr< BoundaryEle > boundary_ele_rhs
MoFEMErrorCode set_blockData(std::map< int, BlockData > &block_data_map)
Core (interface) class.
Definition: Core.hpp:50
MoFEMErrorCode update_slow_rhs(std::string mass_fiedl, boost::shared_ptr< VectorDouble > &mass_ptr)
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1053
EdgeElementForcesAndSourcesCoreBase BoundaryEle
MoFEMErrorCode push_slow_rhs()
PetscErrorCode DMMoFEMTSSetIJacobian(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS Jacobian evaluation function
Definition: DMMMoFEM.cpp:772
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVectorsets the routine to create a global vector associated with the shell DM.
Definition: DMMMoFEM.cpp:1026
MoFEMErrorCode push_stiff_lhs()
MoFEMErrorCode update_vol_fe(boost::shared_ptr< FaceEle > &vol_ele, boost::shared_ptr< PreviousData > &data)
double operator()(const double x, const double y, const double t) const
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
MoFEMErrorCode post_proc_fields()
PetscErrorCode DMMoFEMTSSetRHSFunction(DM dm, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
set TS the right hand side function
Definition: DMMMoFEM.cpp:801
std::vector< boost::shared_ptr< VectorDouble > > dots_ptr
boost::shared_ptr< Ele > vol_ele_slow_rhs
double operator()(const double x, const double y, const double t) const
boost::shared_ptr< Ele > vol_mass_ele
MoFEMErrorCode push_stiff_rhs()
Calculate inverse of jacobian for face element.
Basic algebra on fields.
Definition: FieldBlas.hpp:34
moab::Interface & moab
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
SmartPetscObj< Mat > mass_matrix
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
MoFEMErrorCode set_integration_rule()
MoFEMErrorCode run_analysis()
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
#define CHKERR
Inline error check.
Definition: definitions.h:602
MoFEMErrorCode update_stiff_rhs()
std::vector< Range > inner_surface
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMMoFEM.cpp:915
const double ramp_t
Transform local reference derivatives of shape functions to global derivatives.
constexpr int order
MoFEMErrorCode output_result()
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:291
auto createKSP
Definition: AuxPETSc.hpp:289
Postprocess on face.
MoFEMErrorCode solve()
MoFEMErrorCode set_initial_values(std::string field_name, int block_id, Range &surface, double &init_val)
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
PetscErrorCode DMMoFEMTSSetIFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set TS implicit function evaluation function
Definition: DMMMoFEM.cpp:719
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
continuous field
Definition: definitions.h:177
const double sml
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
Get value at integration points for scalar field.
FTensor::Tensor1< double, 3 > operator()(const double x, const double y, const double t) const
OpCalculateScalarFieldValuesDot OpCalculateScalarValuesDot
boost::shared_ptr< Ele > vol_ele_stiff_rhs
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:67
SmartPetscObj< KSP > mass_Ksp
std::vector< boost::shared_ptr< VectorDouble > > values_ptr
FTensor::Index< 'i', 2 > i
[Common data]
Definition: ContactOps.hpp:26
static char help[]
field with C-1 continuity
Definition: definitions.h:180