v0.10.0
mixed_reac_diff.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
3 #include <RDOperators.hpp>
5 
6 using namespace MoFEM;
7 using namespace ReactionDiffusion;
8 using namespace ErrorEstimates;
9 
10 static char help[] = "...\n\n";
11 
12 // #define M_PI 3.14159265358979323846 /* pi */
13 
14 struct RDProblem {
15 public:
16  RDProblem(MoFEM::Core &core, const int order)
17  : m_field(core)
18  , order(order)
19  , cOmm(m_field.get_comm())
20  , rAnk(m_field.get_comm_rank()) {
21  vol_ele_slow_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
22  natural_bdry_ele_slow_rhs =
23  boost::shared_ptr<BoundaryEle>(new BoundaryEle(m_field));
24  vol_ele_stiff_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
25  vol_ele_stiff_lhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
26 
27  skeleton_fe = boost::shared_ptr<EdgeEle>(new EdgeEle(m_field));
28 
29  post_proc = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
30  new PostProcFaceOnRefinedMesh(m_field));
31 
32  data1 = boost::shared_ptr<PreviousData>(new PreviousData());
33  data2 = boost::shared_ptr<PreviousData>(new PreviousData());
34  data3 = boost::shared_ptr<PreviousData>(new PreviousData());
35 
36  flux_values_ptr1 =
37  boost::shared_ptr<MatrixDouble>(data1, &data1->flux_values);
38 
39  flux_divs_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->flux_divs);
40 
41  mass_values_ptr1 =
42  boost::shared_ptr<VectorDouble>(data1, &data1->mass_values);
43 
44  mass_dots_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->mass_dots);
45 
46  flux_values_ptr2 =
47  boost::shared_ptr<MatrixDouble>(data2, &data2->flux_values);
48  flux_divs_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->flux_divs);
49 
50  mass_values_ptr2 =
51  boost::shared_ptr<VectorDouble>(data2, &data2->mass_values);
52  mass_dots_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->mass_dots);
53 
54  flux_values_ptr3 =
55  boost::shared_ptr<MatrixDouble>(data3, &data3->flux_values);
56  flux_divs_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->flux_divs);
57 
58  mass_values_ptr3 =
59  boost::shared_ptr<VectorDouble>(data3, &data3->mass_values);
60 
61  mass_dots_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->mass_dots);
62  }
63 
64  // RDProblem(const int order) : order(order){}
65  MoFEMErrorCode run_analysis(int nb_sp);
66 
67  double global_error;
68 
69 private:
70  MoFEMErrorCode setup_system();
71  MoFEMErrorCode add_fe(std::string mass_field, std::string flux_field);
72  MoFEMErrorCode set_blockData(std::map<int, BlockData> &block_data_map);
73  MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL, std::string internal);
74  MoFEMErrorCode extract_initial_ents(int block_id, Range &surface);
75  MoFEMErrorCode update_slow_rhs(std::string mass_fiedl,
76  boost::shared_ptr<VectorDouble> &mass_ptr);
77  MoFEMErrorCode push_slow_rhs(std::string mass_field, std::string flux_field,
78  boost::shared_ptr<PreviousData> &data);
79  MoFEMErrorCode update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
80  boost::shared_ptr<PreviousData> &data);
82  update_stiff_rhs(std::string mass_field, std::string flux_field,
83  boost::shared_ptr<VectorDouble> &mass_ptr,
84  boost::shared_ptr<MatrixDouble> &flux_ptr,
85  boost::shared_ptr<VectorDouble> &mass_dot_ptr,
86  boost::shared_ptr<VectorDouble> &flux_div_ptr);
87  MoFEMErrorCode push_stiff_rhs(std::string mass_field, std::string flux_field,
88  boost::shared_ptr<PreviousData> &data,
89  std::map<int, BlockData> &block_map);
90  MoFEMErrorCode update_stiff_lhs(std::string mass_fiedl,
91  std::string flux_field,
92  boost::shared_ptr<VectorDouble> &mass_ptr,
93  boost::shared_ptr<MatrixDouble> &flux_ptr);
94  MoFEMErrorCode push_stiff_lhs(std::string mass_field, std::string flux_field,
95  boost::shared_ptr<PreviousData> &data,
96  std::map<int, BlockData> &block_map);
97 
98  MoFEMErrorCode set_integration_rule();
99  MoFEMErrorCode apply_IC(std::string mass_field, Range &surface,
100  boost::shared_ptr<FaceEle> &initial_ele);
101  MoFEMErrorCode apply_BC(std::string flux_field);
102  MoFEMErrorCode loop_fe();
103  MoFEMErrorCode post_proc_fields(std::string mass_field,
104  std::string flux_field);
105  MoFEMErrorCode output_result();
106  MoFEMErrorCode solve();
107 
108  MoFEM::Interface &m_field;
109  Simple *simple_interface;
112 
113  Range essential_bdry_ents;
114  Range natural_bdry_ents;
115 
117 
118  Range inner_surface1; // nb_species times
119  Range inner_surface2;
120  Range inner_surface3;
121 
122  MPI_Comm cOmm;
123  const int rAnk;
124 
125  int order;
126  int nb_species;
127 
128  std::map<int, BlockData> material_blocks;
129 
130  boost::shared_ptr<FaceEle> vol_ele_slow_rhs;
131  boost::shared_ptr<FaceEle> vol_ele_stiff_rhs;
132  boost::shared_ptr<FaceEle> vol_ele_stiff_lhs;
133  boost::shared_ptr<BoundaryEle> natural_bdry_ele_slow_rhs;
134 
135  boost::shared_ptr<EdgeEle> skeleton_fe;
136 
137  boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc;
138  boost::shared_ptr<Monitor> monitor_ptr;
139 
140  boost::shared_ptr<PreviousData> data1; // nb_species times
141  boost::shared_ptr<PreviousData> data2;
142  boost::shared_ptr<PreviousData> data3;
143 
144  boost::shared_ptr<MatrixDouble> flux_values_ptr1; // nb_species times
145  boost::shared_ptr<MatrixDouble> flux_values_ptr2;
146  boost::shared_ptr<MatrixDouble> flux_values_ptr3;
147 
148  boost::shared_ptr<VectorDouble> flux_divs_ptr1; // nb_species times
149  boost::shared_ptr<VectorDouble> flux_divs_ptr2;
150  boost::shared_ptr<VectorDouble> flux_divs_ptr3;
151 
152  boost::shared_ptr<VectorDouble> mass_values_ptr1; // nb_species times
153  boost::shared_ptr<VectorDouble> mass_values_ptr2;
154  boost::shared_ptr<VectorDouble> mass_values_ptr3;
155 
156  boost::shared_ptr<VectorDouble> mass_dots_ptr1; // nb_species times
157  boost::shared_ptr<VectorDouble> mass_dots_ptr2;
158  boost::shared_ptr<VectorDouble> mass_dots_ptr3;
159 
160  boost::shared_ptr<ForcesAndSourcesCore> null;
161 };
162 
163 const double ramp_t = 1.0;
164 const double sml = 0.0;
165 const double T = M_PI / 2.0;
166 const double R = 0.75;
167 
168 struct KinkFunction {
169  double operator() (const double x) const {
170  return 1 - abs(x);
171  }
172 };
173 
175  double operator()(const double x) const {
176  if(x > 0)
177  {return -1;}
178  else
179  {return 1;}
180 }
181 };
182 
184  double operator()(const double x, const double y, const double t) const {
185  // double g = cos(T * x) * cos(T * y);
186  double g = 0;
187  double r = x * x + y * y;
188  double d = R * R - r;
189 
190  if(d > 0.0){
191  g = exp(-r / d);
192  }
193  if (t <= ramp_t) {
194  return g * t;
195  } else {
196  return g * ramp_t;
197  }
198  }
199 };
200 
202  FTensor::Tensor1<double, 3> operator()(const double x, const double y,
203  const double t) const {
205  // double mx = - T * sin(T * x) * cos(T * y);
206  // double my = - T * cos(T * x) * sin(T * y);
207 
208  double mx = 0.0;
209  double my = 0.0;
210  double r = x * x + y * y;
211  double d = R * R - r;
212 
213  if (d > 0.0) {
214  double rx = 2.0 * x;
215  double ry = 2.0 * y;
216 
217  double g = exp(-r / d);
218 
219  mx = g * (- rx * R * R / (d * d));
220  my = g * (- ry * R * R / (d * d));
221  }
222 
223  if (t <= ramp_t) {
224  grad(0) = mx * t;
225  grad(1) = my * t;
226  } else {
227  grad(0) = mx * ramp_t;
228  grad(1) = my * ramp_t;
229  }
230  grad(2) = 0.0;
231  return grad;
232  }
233 };
234 
236  double operator()(const double x, const double y, const double t) const {
237  // double glap = -2.0 * pow(T, 2) * cos(T * x) * cos(T * y);
238  double glap = 0.0;
239  double r = x * x + y * y;
240  double d = R * R - r;
241  if (d > 0.0) {
242  double rx = 2.0 * x;
243  double ry = 2.0 * y;
244 
245  double rxx = 2.0;
246  double ryy = 2.0;
247  double g = exp(-r / d);
248  double cof = R * R / (d * d);
249  glap = g * ((-rx * cof) * (-rx * cof) + (-ry * cof) * (-ry * cof)
250  - R * R * (rxx * d + 2.0 * rx * rx) / (d * d * d)
251  - R * R * (ryy * d + 2.0 * ry * ry) / (d * d * d));
252  }
253 
254  if (t <= ramp_t) {
255  return glap * t;
256  } else {
257  return glap * ramp_t;
258  }
259  }
260 };
261 
263  double operator()(const double x, const double y, const double t) const {
264  // double gdot = cos(T * x) * cos(T * y);
265 
266  double gdot = 0;
267  double r = x * x + y * y;
268  double d = R * R - r;
269 
270  if (d > 0.0) {
271  gdot = exp(-r / d);
272  }
273  if (t <= ramp_t) {
274  return gdot;
275  } else {
276  return 0;
277  }
278  }
279 };
280 
283  CHKERR m_field.getInterface(simple_interface);
284  CHKERR simple_interface->getOptions();
285  CHKERR simple_interface->loadFile();
287 }
288 
289 MoFEMErrorCode RDProblem::add_fe(std::string mass_field,
290  std::string flux_field) {
292  CHKERR simple_interface->addDomainField(mass_field, L2,
294 
295  CHKERR simple_interface->addDomainField(flux_field, HCURL,
297 
298  CHKERR simple_interface->addBoundaryField(flux_field, HCURL,
300  CHKERR simple_interface->addSkeletonField(flux_field, HCURL,
302 
303  // CHKERR simple_interface->addDataField("ERROR", L2, AINSWORTH_LEGENDRE_BASE,
304  // 1);
305  // CHKERR simple_interface->addDataField("ERROR2", L2, AINSWORTH_LEGENDRE_BASE,
306  // 1);
307 
308  CHKERR simple_interface->setFieldOrder(mass_field, order - 1);
309  CHKERR simple_interface->setFieldOrder(flux_field, order);
310  // CHKERR simple_interface->setFieldOrder("ERROR", 0); // approximation order for error
311  // CHKERR simple_interface->setFieldOrder("ERROR2", 0);
312 
314 }
315 MoFEMErrorCode RDProblem::set_blockData(std::map<int, BlockData> &block_map) {
318  string name = it->getName();
319  const int id = it->getMeshsetId();
320  if (name.compare(0, 14, "REGION1") == 0) {
321  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
322  id, BLOCKSET, 2, block_map[id].block_ents, true);
323  block_map[id].B0 = 1e0;
324  block_map[id].block_id = id;
325  } else if (name.compare(0, 14, "REGION2") == 0) {
326  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
327  id, BLOCKSET, 2, block_map[id].block_ents, true);
328  block_map[id].B0 = 1e0;
329  block_map[id].block_id = id;
330  } else if (name.compare(0, 14, "REGION3") == 0) {
331  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
332  id, BLOCKSET, 2, block_map[id].block_ents, true);
333  block_map[id].B0 = 1e-4;
334  block_map[id].block_id = id;
335  } else if (name.compare(0, 14, "REGION4") == 0) {
336  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
337  id, BLOCKSET, 2, block_map[id].block_ents, true);
338  block_map[id].B0 = 1e-4;
339  block_map[id].block_id = id;
340  } else if (name.compare(0, 14, "REGION5") == 0) {
341  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
342  id, BLOCKSET, 2, block_map[id].block_ents, true);
343  block_map[id].B0 = 1e-4;
344  block_map[id].block_id = id;
345  }
346  }
348 }
349 
351  std::string natural,
352  std::string internal) {
355  string name = it->getName();
356  if (name.compare(0, 14, natural) == 0) {
357 
358  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
359  natural_bdry_ents, true);
360  } else if (name.compare(0, 14, essential) == 0) {
361  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
362  essential_bdry_ents, true);
363  } else if (name.compare(0, 14, internal) == 0) {
364  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
365  internal_edge_ents, true);
366  }
367  }
369 }
370 
371 MoFEMErrorCode RDProblem::extract_initial_ents(int block_id, Range &surface) {
373  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(block_id,
374  BLOCKSET)) {
375  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
376  block_id, BLOCKSET, 2, surface, true);
377  }
379 }
381 RDProblem::update_slow_rhs(std::string mass_field,
382  boost::shared_ptr<VectorDouble> &mass_ptr) {
384  vol_ele_slow_rhs->getOpPtrVector().push_back(
385  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
387 }
388 
390  std::string flux_field,
391  boost::shared_ptr<PreviousData> &data) {
393 
394  // vol_ele_slow_rhs->getOpPtrVector().push_back(
395  // new OpAssembleSlowRhsV(mass_field, data, ExactFunction(), ExactFunctionDot(), ExactFunctionLap()));
396 
397  // natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
398  // new OpSkeletonSource(mass_field,
399  // KinkFunction(),
400  // ExactFunction(),
401  // internal_edge_ents));
402 
403  // natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
404  // new OpAssembleNaturalBCRhsTau(flux_field, natural_bdry_ents));
405 
406  // natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
407  // new OpEssentialBC(flux_field, essential_bdry_ents));
409 }
410 
411 MoFEMErrorCode RDProblem::update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
412  boost::shared_ptr<PreviousData> &data) {
414  vol_ele->getOpPtrVector().push_back(new OpCalculateJacForFace(data->jac));
415  vol_ele->getOpPtrVector().push_back(
416  new OpCalculateInvJacForFace(data->inv_jac));
417  vol_ele->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
418 
419  vol_ele->getOpPtrVector().push_back(
420  new OpSetContravariantPiolaTransformFace(data->jac));
421 
422  vol_ele->getOpPtrVector().push_back(new OpSetInvJacHcurlFace(data->inv_jac));
424 }
426 RDProblem::update_stiff_rhs(std::string mass_field, std::string flux_field,
427  boost::shared_ptr<VectorDouble> &mass_ptr,
428  boost::shared_ptr<MatrixDouble> &flux_ptr,
429  boost::shared_ptr<VectorDouble> &mass_dot_ptr,
430  boost::shared_ptr<VectorDouble> &flux_div_ptr) {
431 
433 
434  vol_ele_stiff_rhs->getOpPtrVector().push_back(
435  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
436 
437  vol_ele_stiff_rhs->getOpPtrVector().push_back(
438  new OpCalculateHdivVectorField<3>(flux_field, flux_ptr));
439 
440  vol_ele_stiff_rhs->getOpPtrVector().push_back(
441  new OpCalculateScalarValuesDot(mass_field, mass_dot_ptr));
442 
443  vol_ele_stiff_rhs->getOpPtrVector().push_back(
444  new OpCalculateHdivVectorDivergence<3, 2>(flux_field, flux_div_ptr));
446 }
447 
449  std::string flux_field,
450  boost::shared_ptr<PreviousData> &data,
451  std::map<int, BlockData> &block_map) {
453  vol_ele_stiff_rhs->getOpPtrVector().push_back(
454  new OpAssembleStiffRhsTau<3>(flux_field, data, block_map));
455 
456  vol_ele_stiff_rhs->getOpPtrVector().push_back(
457  new OpAssembleStiffRhsV<3>(mass_field, data, block_map, ExactFunction(),
460 }
461 
463 RDProblem::update_stiff_lhs(std::string mass_field, std::string flux_field,
464  boost::shared_ptr<VectorDouble> &mass_ptr,
465  boost::shared_ptr<MatrixDouble> &flux_ptr) {
467  vol_ele_stiff_lhs->getOpPtrVector().push_back(
468  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
469 
470  vol_ele_stiff_lhs->getOpPtrVector().push_back(
471  new OpCalculateHdivVectorField<3>(flux_field, flux_ptr));
473 }
474 
476  std::string flux_field,
477  boost::shared_ptr<PreviousData> &data,
478  std::map<int, BlockData> &block_map) {
480  vol_ele_stiff_lhs->getOpPtrVector().push_back(
481  new OpAssembleLhsTauTau<3>(flux_field, data, block_map));
482 
483  vol_ele_stiff_lhs->getOpPtrVector().push_back(
484  new OpAssembleLhsVV(mass_field));
485 
486  vol_ele_stiff_lhs->getOpPtrVector().push_back(
487  new OpAssembleLhsTauV<3>(flux_field, mass_field, data, block_map));
488 
489  vol_ele_stiff_lhs->getOpPtrVector().push_back(
490  new OpAssembleLhsVTau(mass_field, flux_field));
492 }
493 
496  auto vol_rule = [](int, int, int p) -> int {
497  return 2 * p; };
498  vol_ele_slow_rhs->getRuleHook = vol_rule;
499  natural_bdry_ele_slow_rhs->getRuleHook = vol_rule;
500 
501  vol_ele_stiff_rhs->getRuleHook = vol_rule;
502 
503  vol_ele_stiff_lhs->getRuleHook = vol_rule;
504  skeleton_fe->getRuleHook = vol_rule;
506 }
507 
508 MoFEMErrorCode RDProblem::apply_IC(std::string mass_field, Range &surface,
509  boost::shared_ptr<FaceEle> &initial_ele) {
511  initial_ele->getOpPtrVector().push_back(
512  new OpInitialMass(mass_field, surface));
514 }
515 
516 MoFEMErrorCode RDProblem::apply_BC(std::string flux_field) {
518  CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
519  "SimpleProblem", flux_field, essential_bdry_ents);
520 
522 }
525 
526  // Vec R, U;
527  // CHKERR DMCreateGlobalVector(dm, &R);
528 
529  // CHKERR VecDuplicate(R, &U);
530 
531  // DM dm_sub_ff, dm_sub_fm;
532  // ublas::matrix<Mat> nested_matrices(2, 2);
533  // ublas::vector<IS> nested_is(2);
534 
535  // CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_ff);
536  // CHKERR DMSetType(dm_sub_ff, "DMMOFEM");
537  // CHKERR DMMoFEMCreateSubDM(dm_sub_hh, dm, "SUB_FF");
538 
539  // CHKERR DMMoFEMSetSquareProblem(dm_sub_ff, PETSC_TRUE);
540  // CHKERR DMMoFEMAddSubFieldRow(dm_sub_ff, "FLUX1");
541  // CHKERR DMMoFEMAddSubFieldCol(dm_sub_ff, "FLUX1");
542  // CHKERR DMMoFEMAddElement(dm_sub_ff,
543  // simple_interface->getDomainFEName().c_str());
544  // CHKERR DMSetUp(dm_sub_ff);
545 
546  // CHKERR DMMoFEMGetSubRowIS(dm_sub_ff, &nested_is[0]);
547  // CHKERR DMCreateMatrix(dm_sub_ff, &nested_matrices(0, 0));
548 
549  // vol_ele_stiff_lhs->getFEMethod()->ts_B = nested_matrices(0, 0);
550  // CHKERR TSSetType(ts, TSARKIMEX);
551  // CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
552 
553  // CHKERR DMMoFEMTSSetIJacobian(dm_sub_ff, simple_interface->getDomainFEName(),
554  // vol_ele_stiff_lhs, null, null);
555 
556  // CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
557  // CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
558 
559  CHKERR TSSetType(ts, TSARKIMEX);
560  CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
561 
562  CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
563  vol_ele_stiff_lhs, null, null);
564 
565 
566  CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
567  vol_ele_stiff_rhs, null, null);
568 
569  CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
570  vol_ele_slow_rhs, null, null);
571 
572  // CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getSkeletonFEName(),
573  // skeleton_fe, null, null);
574 
576 }
577 
579  std::string flux_field) {
581  post_proc->addFieldValuesPostProc(mass_field);
582  post_proc->addFieldValuesPostProc(flux_field);
584 }
585 
588  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
589  monitor_ptr, null, null);
591 }
594  // Create solution vector
597  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
598  // Solve problem
599  double ftime = 1;
600  CHKERR TSSetDM(ts, dm);
601  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
602  CHKERR TSSetSolution(ts, X);
603  CHKERR TSSetFromOptions(ts);
604 
605  if (0) {
606  SNES snes;
607  CHKERR TSGetSNES(ts, &snes);
608  KSP ksp;
609  CHKERR SNESGetKSP(snes, &ksp);
610  PC pc;
611  CHKERR KSPGetPC(ksp, &pc);
612  PetscBool is_pcfs = PETSC_FALSE;
613  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
614  // Set up FIELDSPLIT
615  // Only is user set -pc_type fieldsplit
616  if (is_pcfs == PETSC_TRUE) {
617  IS is_mass, is_flux;
618  const MoFEM::Problem *problem_ptr;
619  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
620  CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
621  problem_ptr->getName(), ROW, "MASS1", 0, 1, &is_mass);
622  CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
623  problem_ptr->getName(), ROW, "FLUX1", 0, 1, &is_flux);
624  // CHKERR ISView(is_flux, PETSC_VIEWER_STDOUT_SELF);
625  // CHKERR ISView(is_mass, PETSC_VIEWER_STDOUT_SELF);
626 
627  CHKERR PCFieldSplitSetIS(pc, NULL, is_flux);
628  CHKERR PCFieldSplitSetIS(pc, NULL, is_mass);
629 
630 
631  CHKERR ISDestroy(&is_flux);
632  CHKERR ISDestroy(&is_mass);
633  }
634  }
635 
636  CHKERR TSSolve(ts, X);
638 }
639 
642 
643  // set nb_species
644  CHKERR setup_system(); // only once
645  nb_species = nb_sp;
646  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
647  CHKERR add_fe("MASS1", "FLUX1"); // nb_species times
648  if (nb_species == 2 || nb_species == 3) {
649  CHKERR add_fe("MASS2", "FLUX2");
650  if (nb_species == 3) {
651  CHKERR add_fe("MASS3", "FLUX3");
652  }
653  }
654  }
655 
656  CHKERR simple_interface->setUp();
657 
658  CHKERR set_blockData(material_blocks);
659 
660  CHKERR extract_bd_ents("ESSENTIAL", "NATURAL", "INTERNAL"); // nb_species times
661 
662  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
663  CHKERR extract_initial_ents(2, inner_surface1);
664  CHKERR update_slow_rhs("MASS1", mass_values_ptr1);
665  if (nb_species == 1) {
666  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
667  "MASS1", data1, data1, data1, material_blocks));
668  } else if (nb_species == 2 || nb_species == 3) {
669  CHKERR extract_initial_ents(3, inner_surface2);
670  CHKERR update_slow_rhs("MASS2", mass_values_ptr2);
671  if (nb_species == 2) {
672  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
673  "MASS1", data1, data2, data2, material_blocks));
674  } else if (nb_species == 3) {
675  CHKERR extract_initial_ents(4, inner_surface3);
676  CHKERR update_slow_rhs("MASS3", mass_values_ptr3);
677  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
678  "MASS1", data1, data2, data3, material_blocks));
679  }
680  }
681  }
682  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
684 
685  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
686  CHKERR push_slow_rhs("MASS1", "FLUX1", data1); // nb_species times
687  if (nb_species == 2 || nb_species == 3) {
688  CHKERR push_slow_rhs("MASS2", "FLUX2", data2);
689  if (nb_species == 3) {
690  CHKERR push_slow_rhs("MASS3", "FLUX3", data3);
691  }
692  }
693  }
694 
695  CHKERR update_vol_fe(vol_ele_stiff_rhs, data1);
696 
697  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
698  CHKERR update_stiff_rhs("MASS1", "FLUX1", mass_values_ptr1,
699  flux_values_ptr1, mass_dots_ptr1, flux_divs_ptr1);
700  CHKERR push_stiff_rhs("MASS1", "FLUX1", data1,
701  material_blocks); // nb_species times
702  if (nb_species == 2 || nb_species == 3) {
703  CHKERR update_stiff_rhs("MASS2", "FLUX2", mass_values_ptr2,
704  flux_values_ptr2, mass_dots_ptr2, flux_divs_ptr2);
705  CHKERR push_stiff_rhs("MASS2", "FLUX2", data2, material_blocks);
706  if (nb_species == 3) {
707  CHKERR update_stiff_rhs("MASS3", "FLUX3", mass_values_ptr3,
708  flux_values_ptr3, mass_dots_ptr3,
709  flux_divs_ptr3);
710  CHKERR push_stiff_rhs("MASS3", "FLUX3", data3, material_blocks);
711  }
712  }
713  }
714 
715  CHKERR update_vol_fe(vol_ele_stiff_lhs, data1);
716 
717  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
718  CHKERR update_stiff_lhs("MASS1", "FLUX1", mass_values_ptr1,
719  flux_values_ptr1);
720  CHKERR push_stiff_lhs("MASS1", "FLUX1", data1,
721  material_blocks); // nb_species times
722  if (nb_species == 2 || nb_species == 3) {
723  CHKERR update_stiff_lhs("MASS2", "FLUX2", mass_values_ptr2,
724  flux_values_ptr2);
725  CHKERR push_stiff_lhs("MASS2", "FLUX2", data2, material_blocks);
726  if (nb_species == 3) {
727  CHKERR update_stiff_lhs("MASS3", "FLUX3", mass_values_ptr3,
728  flux_values_ptr3);
729  CHKERR push_stiff_lhs("MASS3", "FLUX3", data3, material_blocks);
730  }
731  }
732  }
733  global_error = 0;
734  // vol_ele_stiff_rhs->getOpPtrVector().push_back(
735  // new OpError(ExactFunction(), ExactFunctionLap(), ExactFunctionGrad(),
736  // data1, material_blocks, global_error));
737 
738  JumpData jump_data;
739 
740  // boost::shared_ptr<EdgeEle> skeleton_fe =
741  // boost::shared_ptr<EdgeEle>(new EdgeEle(m_field));
742 
743  // skeleton_fe->getOpPtrVector().push_back(
744  // new OpSetContrariantPiolaTransformOnEdge());
745 
746  skeleton_fe->getOpPtrVector().push_back(new OpCalculateJumpErrors(
747  "FLUX1", "MASS1", "ERROR2", m_field, jump_data));
748  CHKERR set_integration_rule();
749  dm = simple_interface->getDM();
750  ts = createTS(m_field.get_comm());
751  // boost::shared_ptr<FaceEle> initial_mass_ele(new FaceEle(m_field));
752 
753  // if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
754  // CHKERR apply_IC("MASS1", inner_surface1,
755  // initial_mass_ele); // nb_species times
756  // if (nb_species == 2 || nb_species == 3) {
757  // CHKERR apply_IC("MASS2", inner_surface2, initial_mass_ele);
758  // if (nb_species == 3) {
759  // CHKERR apply_IC("MASS3", inner_surface3, initial_mass_ele);
760  // }
761  // }
762  // }
763  // CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
764  // initial_mass_ele);
765 
766  // if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
767  // CHKERR apply_BC("FLUX1"); // nb_species times
768  // if (nb_species == 2 || nb_species == 3) {
769  // CHKERR apply_BC("FLUX2");
770  // if (nb_species == 3) {
771  // CHKERR apply_BC("FLUX3");
772  // }
773  // }
774  // }
775 
776  CHKERR loop_fe(); // only once
777  post_proc->generateReferenceElementMesh(); // only once
778 
779 
780 
781  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
782  CHKERR post_proc_fields("MASS1", "FLUX1");
783  // post_proc->addFieldValuesPostProc("ERROR");
784  if (nb_species == 2 || nb_species == 3) {
785  CHKERR post_proc_fields("MASS2", "FLUX2");
786  if (nb_species == 3) {
787  CHKERR post_proc_fields("MASS3", "FLUX3");
788  }
789  }
790  }
791 
792  //the double global_vector is a global_vector which sums errors in each element
793 
794  // Vec gVec;
795  // CHKERR VecCreateMPI(PETSC_COMM_WORLD, 1, 1, &gVec);
796  // double err_per_step = 0;
797 
798  // auto compute_global_error = [&gVec](double &g_err, double &err_out) {
799  // MoFEMFunctionBegin;
800 
801  // CHKERR VecZeroEntries(gVec);
802  // CHKERR VecAssemblyBegin(gVec);
803  // CHKERR VecAssemblyEnd(gVec);
804 
805  // int ind[1] = {0};
806  // CHKERR VecSetValues(gVec, 1, ind, &err_out, ADD_VALUES);
807  // CHKERR VecAssemblyBegin(gVec);
808  // CHKERR VecAssemblyEnd(gVec);
809 
810  // CHKERR VecGetValues(gVec, 1, ind, &err_out);
811 
812  // MoFEMFunctionReturn(0);
813  // };
814 
815  // Vec error_per_proc;
816  // CHKERR VecCreateMPI(cOmm, 1, PETSC_DECIDE, &error_per_proc);
817  // auto get_global_error = [&]() {
818  // MoFEMFunctionBegin;
819  // CHKERR VecSetValue(error_per_proc, m_field.get_comm_rank(), global_error, INSERT_VALUES);
820  // MoFEMFunctionReturn(0);
821  // };
822 
823  // auto reinitialize_global = [&]() {
824  // MoFEMFunctionBegin;
825  // CHKERR get_global_error();
826  // CHKERR VecAssemblyBegin(error_per_proc);
827  // CHKERR VecAssemblyEnd(error_per_proc);
828  // double error_sum;
829  // CHKERR VecSum(error_per_proc, &error_sum);
830  // CHKERR PetscPrintf(PETSC_COMM_WORLD, "Error : %3.4e \n",
831  // error_sum);
832  // global_error = 0;
833  // MoFEMFunctionReturn(0);
834  // };
835 
836  // vol_ele_stiff_lhs->postProcessHook = reinitialize_global;
837 
838 
839 
840  monitor_ptr = boost::shared_ptr<Monitor>(
841  new Monitor(cOmm, rAnk, dm, post_proc, global_error)); // nb_species times
842  CHKERR output_result(); // only once
843  CHKERR solve(); // only once
845 }
846 
847 int main(int argc, char *argv[]) {
848  const char param_file[] = "param_file.petsc";
849  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
850  try {
851  moab::Core mb_instance;
852  moab::Interface &moab = mb_instance;
853  MoFEM::Core core(moab);
854  DMType dm_name = "DMMOFEM";
855  CHKERR DMRegister_MoFEM(dm_name);
856 
857  int order = 1;
858  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
859  int nb_species = 1;
860  RDProblem reac_diff_problem(core, order);
861  CHKERR reac_diff_problem.run_analysis(nb_species);
862  }
863  CATCH_ERRORS;
865  return 0;
866 }
Calculate divergence of vector field.
const double sml
boost::shared_ptr< VectorDouble > flux_divs_ptr3
Deprecated interface functions.
Problem manager is used to build and partition problems.
MoFEMErrorCode add_fe()
const double ramp_t
MoFEMErrorCode update_stiff_lhs()
CoreTmp< 0 > Core
Definition: Core.hpp:1129
MoFEMErrorCode setup_system()
auto createTS
Definition: AuxPETSc.hpp:275
Section manager is used to create indexes and sectionsFIXME: ISManager is not properly testsed by ato...
Definition: ISManager.hpp:36
boost::shared_ptr< MatrixDouble > flux_values_ptr3
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
auto post_proc
double operator()(const double x, const double y, const double t) const
static Index< 'p', 3 > p
MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
int main(int argc, char *argv[])
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode set_blockData(std::map< int, BlockData > &block_data_map)
MoFEMErrorCode update_slow_rhs(std::string mass_fiedl, boost::shared_ptr< VectorDouble > &mass_ptr)
keeps basic data about problemThis is low level structure with information about problem,...
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
MoFEMErrorCode apply_BC(std::string flux_field)
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
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
double operator()(const double x, const double y, const double t) const
MoFEMErrorCode push_stiff_rhs()
RDProblem(MoFEM::Core &core, const int order)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
MoFEM::EdgeElementForcesAndSourcesCoreSwitch< EdgeElementForcesAndSourcesCore::NO_HO_GEOMETRY > EdgeEle
constexpr int order
MoFEMErrorCode extract_initial_ents(int block_id, Range &surface)
Calculate inverse of jacobian for face element.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Exact gradient.
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
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()
boost::shared_ptr< EdgeEle > skeleton_fe
MoFEMErrorCode loop_fe()
field with continuous tangents
Definition: definitions.h:178
Get vector field for H-div approximation.
boost::shared_ptr< VectorDouble > mass_values_ptr3
MoFEMErrorCode run_analysis()
boost::shared_ptr< PreviousData > data3
#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:604
MoFEMErrorCode update_stiff_rhs()
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 r
rate factor
Range internal_edge_ents
const double R
MoFEMErrorCode output_result()
Postprocess on face.
MoFEMErrorCode solve()
brief Transform local reference derivatives of shape function to global derivatives for face
MoFEMErrorCode apply_IC(std::string mass_field, Range &surface, boost::shared_ptr< FaceEle > &initial_ele, double &init_val)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
Apply contravariant (Piola) transfer to Hdiv space on face.
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
Core (interface) class.
Definition: Core.hpp:77
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
Calculate jacobian for face element.
std::string getName() const
std::string param_file
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
double operator()(const double x) const
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
boost::shared_ptr< VectorDouble > mass_dots_ptr3
static char help[]
Make Hdiv space from Hcurl space in 2d.
field with C-1 continuity
Definition: definitions.h:180