v0.10.0
mixed_reac_diff.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <BasicFiniteElements.hpp>
3 #include <RDOperators.hpp>
4 
5 using namespace MoFEM;
6 using namespace ReactionDiffusion;
7 
8 static char help[] = "...\n\n";
9 
10 // #define M_PI 3.14159265358979323846 /* pi */
11 
12 struct RDProblem {
13 public:
14  RDProblem(MoFEM::Core &core, const int order)
15  : m_field(core)
16  , order(order)
17  , cOmm(m_field.get_comm())
18  , rAnk(m_field.get_comm_rank()) {
19  vol_ele_slow_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
20  natural_bdry_ele_slow_rhs =
21  boost::shared_ptr<BoundaryEle>(new BoundaryEle(m_field));
22  vol_ele_stiff_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
23  vol_ele_stiff_lhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
24  post_proc = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
25  new PostProcFaceOnRefinedMesh(m_field));
26 
27  data1 = boost::shared_ptr<PreviousData>(new PreviousData());
28  data2 = boost::shared_ptr<PreviousData>(new PreviousData());
29  data3 = boost::shared_ptr<PreviousData>(new PreviousData());
30 
31  flux_values_ptr1 =
32  boost::shared_ptr<MatrixDouble>(data1, &data1->flux_values);
33 
34  flux_divs_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->flux_divs);
35 
36  mass_values_ptr1 =
37  boost::shared_ptr<VectorDouble>(data1, &data1->mass_values);
38 
39  mass_dots_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->mass_dots);
40 
41  flux_values_ptr2 =
42  boost::shared_ptr<MatrixDouble>(data2, &data2->flux_values);
43  flux_divs_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->flux_divs);
44 
45  mass_values_ptr2 =
46  boost::shared_ptr<VectorDouble>(data2, &data2->mass_values);
47  mass_dots_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->mass_dots);
48 
49  flux_values_ptr3 =
50  boost::shared_ptr<MatrixDouble>(data3, &data3->flux_values);
51  flux_divs_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->flux_divs);
52 
53  mass_values_ptr3 =
54  boost::shared_ptr<VectorDouble>(data3, &data3->mass_values);
55 
56  mass_dots_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->mass_dots);
57  }
58 
59  // RDProblem(const int order) : order(order){}
60  MoFEMErrorCode run_analysis(int nb_sp);
61 
62  double global_error;
63 
64 private:
66  MoFEMErrorCode add_fe(std::string mass_field, std::string flux_field);
67  MoFEMErrorCode set_blockData(std::map<int, BlockData> &block_data_map);
68  MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL, std::string internal);
69  MoFEMErrorCode extract_initial_ents(int block_id, Range &surface);
70  MoFEMErrorCode update_slow_rhs(std::string mass_fiedl,
71  boost::shared_ptr<VectorDouble> &mass_ptr);
72  MoFEMErrorCode push_slow_rhs(std::string mass_field, std::string flux_field,
73  boost::shared_ptr<PreviousData> &data);
74  MoFEMErrorCode update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
75  boost::shared_ptr<PreviousData> &data);
77  update_stiff_rhs(std::string mass_field, std::string flux_field,
78  boost::shared_ptr<VectorDouble> &mass_ptr,
79  boost::shared_ptr<MatrixDouble> &flux_ptr,
80  boost::shared_ptr<VectorDouble> &mass_dot_ptr,
81  boost::shared_ptr<VectorDouble> &flux_div_ptr);
82  MoFEMErrorCode push_stiff_rhs(std::string mass_field, std::string flux_field,
83  boost::shared_ptr<PreviousData> &data,
84  std::map<int, BlockData> &block_map);
85  MoFEMErrorCode update_stiff_lhs(std::string mass_fiedl,
86  std::string flux_field,
87  boost::shared_ptr<VectorDouble> &mass_ptr,
88  boost::shared_ptr<MatrixDouble> &flux_ptr);
89  MoFEMErrorCode push_stiff_lhs(std::string mass_field, std::string flux_field,
90  boost::shared_ptr<PreviousData> &data,
91  std::map<int, BlockData> &block_map);
92 
94  MoFEMErrorCode apply_IC(std::string mass_field, Range &surface,
95  boost::shared_ptr<FaceEle> &initial_ele);
96  MoFEMErrorCode apply_BC(std::string flux_field);
98  MoFEMErrorCode post_proc_fields(std::string mass_field,
99  std::string flux_field);
102 
103  MoFEM::Interface &m_field;
104  Simple *simple_interface;
107 
108  Range essential_bdry_ents;
109  Range natural_bdry_ents;
110 
112 
113  Range inner_surface1; // nb_species times
114  Range inner_surface2;
115  Range inner_surface3;
116 
117  MPI_Comm cOmm;
118  const int rAnk;
119 
120  int order;
121  int nb_species;
122 
123  std::map<int, BlockData> material_blocks;
124 
125  boost::shared_ptr<FaceEle> vol_ele_slow_rhs;
126  boost::shared_ptr<FaceEle> vol_ele_stiff_rhs;
127  boost::shared_ptr<FaceEle> vol_ele_stiff_lhs;
128  boost::shared_ptr<BoundaryEle> natural_bdry_ele_slow_rhs;
129 
130  boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc;
131  boost::shared_ptr<Monitor> monitor_ptr;
132 
133  boost::shared_ptr<PreviousData> data1; // nb_species times
134  boost::shared_ptr<PreviousData> data2;
135  boost::shared_ptr<PreviousData> data3;
136 
137  boost::shared_ptr<MatrixDouble> flux_values_ptr1; // nb_species times
138  boost::shared_ptr<MatrixDouble> flux_values_ptr2;
139  boost::shared_ptr<MatrixDouble> flux_values_ptr3;
140 
141  boost::shared_ptr<VectorDouble> flux_divs_ptr1; // nb_species times
142  boost::shared_ptr<VectorDouble> flux_divs_ptr2;
143  boost::shared_ptr<VectorDouble> flux_divs_ptr3;
144 
145  boost::shared_ptr<VectorDouble> mass_values_ptr1; // nb_species times
146  boost::shared_ptr<VectorDouble> mass_values_ptr2;
147  boost::shared_ptr<VectorDouble> mass_values_ptr3;
148 
149  boost::shared_ptr<VectorDouble> mass_dots_ptr1; // nb_species times
150  boost::shared_ptr<VectorDouble> mass_dots_ptr2;
151  boost::shared_ptr<VectorDouble> mass_dots_ptr3;
152 
153  boost::shared_ptr<ForcesAndSourcesCore> null;
154 };
155 
156 const double ramp_t = 1.0;
157 const double sml = 0.0;
158 const double T = M_PI / 2.0;
159 
160 struct KinkFunction {
161  double operator() (const double x) const {
162  return 1 - abs(x);
163  }
164 };
165 
167  double operator()(const double x) const {
168  if(x > 0)
169  {return -1;}
170  else
171  {return 1;}
172 }
173 };
174 
176  double operator()(const double x, const double y, const double t) const {
177  double g = cos(T * x) * cos(T * y);
178  if (t <= ramp_t) {
179  return g * t;
180  } else {
181  return g * ramp_t;
182  }
183  }
184 };
185 
187  FTensor::Tensor1<double, 3> operator()(const double x, const double y,
188  const double t) const {
190  double mx = - T * sin(T * x) * cos(T * y);
191  double my = - T * cos(T * x) * sin(T * y);
192 
193 
194  if (t <= ramp_t) {
195  grad(0) = mx * t;
196  grad(1) = my * t;
197  } else {
198  grad(0) = mx * ramp_t;
199  grad(1) = my * ramp_t;
200  }
201  grad(2) = 0.0;
202  return grad;
203  }
204 };
205 
207  double operator()(const double x, const double y, const double t) const {
208  double glap = -2.0 * pow(T, 2) * cos(T * x) * cos(T * y);
209  if (t <= ramp_t) {
210  return glap * t;
211  } else {
212  return glap * ramp_t;
213  }
214  }
215 };
216 
218  double operator()(const double x, const double y, const double t) const {
219  double gdot = cos(T * x) * cos(T * y);
220  if (t <= ramp_t) {
221  return gdot;
222  } else {
223  return 0;
224  }
225  }
226 };
227 
230  CHKERR m_field.getInterface(simple_interface);
231  CHKERR simple_interface->getOptions();
232  CHKERR simple_interface->loadFile();
234 }
235 
236 MoFEMErrorCode RDProblem::add_fe(std::string mass_field,
237  std::string flux_field) {
239  CHKERR simple_interface->addDomainField(mass_field, L2,
241 
242  CHKERR simple_interface->addDomainField(flux_field, HCURL,
244 
245  CHKERR simple_interface->addBoundaryField(flux_field, HCURL,
247  CHKERR simple_interface->addDataField("ERROR", L2, AINSWORTH_LEGENDRE_BASE,
248  1);
249 
250  CHKERR simple_interface->setFieldOrder(mass_field, order - 1);
251  CHKERR simple_interface->setFieldOrder(flux_field, order);
252  CHKERR simple_interface->setFieldOrder("ERROR", 0); // approximation order for error
253 
255 }
256 MoFEMErrorCode RDProblem::set_blockData(std::map<int, BlockData> &block_map) {
259  string name = it->getName();
260  const int id = it->getMeshsetId();
261  if (name.compare(0, 14, "REGION1") == 0) {
262  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
263  id, BLOCKSET, 2, block_map[id].block_ents, true);
264  block_map[id].B0 = 1e-4;
265  block_map[id].block_id = id;
266  } else if (name.compare(0, 14, "REGION2") == 0) {
267  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
268  id, BLOCKSET, 2, block_map[id].block_ents, true);
269  block_map[id].B0 = 1e-4;
270  block_map[id].block_id = id;
271  } else if (name.compare(0, 14, "REGION3") == 0) {
272  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
273  id, BLOCKSET, 2, block_map[id].block_ents, true);
274  block_map[id].B0 = 1e-4;
275  block_map[id].block_id = id;
276  } else if (name.compare(0, 14, "REGION4") == 0) {
277  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
278  id, BLOCKSET, 2, block_map[id].block_ents, true);
279  block_map[id].B0 = 1e-4;
280  block_map[id].block_id = id;
281  } else if (name.compare(0, 14, "REGION5") == 0) {
282  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
283  id, BLOCKSET, 2, block_map[id].block_ents, true);
284  block_map[id].B0 = 1e-4;
285  block_map[id].block_id = id;
286  }
287  }
289 }
290 
292  std::string natural,
293  std::string internal) {
296  string name = it->getName();
297  if (name.compare(0, 14, natural) == 0) {
298 
299  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
300  natural_bdry_ents, true);
301  } else if (name.compare(0, 14, essential) == 0) {
302  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
303  essential_bdry_ents, true);
304  } else if (name.compare(0, 14, internal) == 0) {
305  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
306  internal_edge_ents, true);
307  }
308  }
310 }
311 
312 MoFEMErrorCode RDProblem::extract_initial_ents(int block_id, Range &surface) {
314  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(block_id,
315  BLOCKSET)) {
316  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
317  block_id, BLOCKSET, 2, surface, true);
318  }
320 }
322 RDProblem::update_slow_rhs(std::string mass_field,
323  boost::shared_ptr<VectorDouble> &mass_ptr) {
325  vol_ele_slow_rhs->getOpPtrVector().push_back(
326  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
328 }
329 
331  std::string flux_field,
332  boost::shared_ptr<PreviousData> &data) {
334 
335  vol_ele_slow_rhs->getOpPtrVector().push_back(
336  new OpAssembleSlowRhsV(mass_field, data, ExactFunction(), ExactFunctionDot(), ExactFunctionLap()));
337 
338  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
339  new OpSkeletonSource(mass_field,
340  KinkFunction(),
341  ExactFunction(),
342  internal_edge_ents));
343 
344  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
345  new OpAssembleNaturalBCRhsTau(flux_field, natural_bdry_ents));
346 
347  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
348  new OpEssentialBC(flux_field, essential_bdry_ents));
350 }
351 
352 MoFEMErrorCode RDProblem::update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
353  boost::shared_ptr<PreviousData> &data) {
355  vol_ele->getOpPtrVector().push_back(new OpCalculateJacForFace(data->jac));
356  vol_ele->getOpPtrVector().push_back(
357  new OpCalculateInvJacForFace(data->inv_jac));
358  vol_ele->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
359 
360  vol_ele->getOpPtrVector().push_back(
361  new OpSetContravariantPiolaTransformFace(data->jac));
362 
363  vol_ele->getOpPtrVector().push_back(new OpSetInvJacHcurlFace(data->inv_jac));
365 }
367 RDProblem::update_stiff_rhs(std::string mass_field, std::string flux_field,
368  boost::shared_ptr<VectorDouble> &mass_ptr,
369  boost::shared_ptr<MatrixDouble> &flux_ptr,
370  boost::shared_ptr<VectorDouble> &mass_dot_ptr,
371  boost::shared_ptr<VectorDouble> &flux_div_ptr) {
372 
374 
375  vol_ele_stiff_rhs->getOpPtrVector().push_back(
376  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
377 
378  vol_ele_stiff_rhs->getOpPtrVector().push_back(
379  new OpCalculateHdivVectorField<3>(flux_field, flux_ptr));
380 
381  vol_ele_stiff_rhs->getOpPtrVector().push_back(
382  new OpCalculateScalarValuesDot(mass_field, mass_dot_ptr));
383 
384  vol_ele_stiff_rhs->getOpPtrVector().push_back(
385  new OpCalculateHdivVectorDivergence<3, 2>(flux_field, flux_div_ptr));
387 }
388 
390  std::string flux_field,
391  boost::shared_ptr<PreviousData> &data,
392  std::map<int, BlockData> &block_map) {
394  vol_ele_stiff_rhs->getOpPtrVector().push_back(
395  new OpAssembleStiffRhsTau<3>(flux_field, data, block_map));
396 
397  vol_ele_stiff_rhs->getOpPtrVector().push_back(
398  new OpAssembleStiffRhsV<3>(mass_field, data, block_map, ExactFunction(),
401 }
402 
404 RDProblem::update_stiff_lhs(std::string mass_field, std::string flux_field,
405  boost::shared_ptr<VectorDouble> &mass_ptr,
406  boost::shared_ptr<MatrixDouble> &flux_ptr) {
408  vol_ele_stiff_lhs->getOpPtrVector().push_back(
409  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
410 
411  vol_ele_stiff_lhs->getOpPtrVector().push_back(
412  new OpCalculateHdivVectorField<3>(flux_field, flux_ptr));
414 }
415 
417  std::string flux_field,
418  boost::shared_ptr<PreviousData> &data,
419  std::map<int, BlockData> &block_map) {
421  vol_ele_stiff_lhs->getOpPtrVector().push_back(
422  new OpAssembleLhsTauTau<3>(flux_field, data, block_map));
423 
424  vol_ele_stiff_lhs->getOpPtrVector().push_back(
425  new OpAssembleLhsVV(mass_field));
426 
427  vol_ele_stiff_lhs->getOpPtrVector().push_back(
428  new OpAssembleLhsTauV<3>(flux_field, mass_field, data, block_map));
429 
430  vol_ele_stiff_lhs->getOpPtrVector().push_back(
431  new OpAssembleLhsVTau(mass_field, flux_field));
433 }
434 
437  auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
438  vol_ele_slow_rhs->getRuleHook = vol_rule;
439  natural_bdry_ele_slow_rhs->getRuleHook = vol_rule;
440 
441  vol_ele_stiff_rhs->getRuleHook = vol_rule;
442 
443  vol_ele_stiff_lhs->getRuleHook = vol_rule;
445 }
446 
447 MoFEMErrorCode RDProblem::apply_IC(std::string mass_field, Range &surface,
448  boost::shared_ptr<FaceEle> &initial_ele) {
450  initial_ele->getOpPtrVector().push_back(
451  new OpInitialMass(mass_field, surface));
453 }
454 
455 MoFEMErrorCode RDProblem::apply_BC(std::string flux_field) {
457  CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
458  "SimpleProblem", flux_field, essential_bdry_ents);
459 
461 }
464 
465  // Vec R, U;
466  // CHKERR DMCreateGlobalVector(dm, &R);
467 
468  // CHKERR VecDuplicate(R, &U);
469 
470  // DM dm_sub_ff, dm_sub_fm;
471  // ublas::matrix<Mat> nested_matrices(2, 2);
472  // ublas::vector<IS> nested_is(2);
473 
474  // CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_ff);
475  // CHKERR DMSetType(dm_sub_ff, "DMMOFEM");
476  // CHKERR DMMoFEMCreateSubDM(dm_sub_hh, dm, "SUB_FF");
477 
478  // CHKERR DMMoFEMSetSquareProblem(dm_sub_ff, PETSC_TRUE);
479  // CHKERR DMMoFEMAddSubFieldRow(dm_sub_ff, "FLUX1");
480  // CHKERR DMMoFEMAddSubFieldCol(dm_sub_ff, "FLUX1");
481  // CHKERR DMMoFEMAddElement(dm_sub_ff,
482  // simple_interface->getDomainFEName().c_str());
483  // CHKERR DMSetUp(dm_sub_ff);
484 
485  // CHKERR DMMoFEMGetSubRowIS(dm_sub_ff, &nested_is[0]);
486  // CHKERR DMCreateMatrix(dm_sub_ff, &nested_matrices(0, 0));
487 
488  // vol_ele_stiff_lhs->getFEMethod()->ts_B = nested_matrices(0, 0);
489  // CHKERR TSSetType(ts, TSARKIMEX);
490  // CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
491 
492  // CHKERR DMMoFEMTSSetIJacobian(dm_sub_ff, simple_interface->getDomainFEName(),
493  // vol_ele_stiff_lhs, null, null);
494 
495  // CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
496  // CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
497 
498  CHKERR TSSetType(ts, TSARKIMEX);
499  CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
500 
501  CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
502  vol_ele_stiff_lhs, null, null);
503 
504 
505  CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
506  vol_ele_stiff_rhs, null, null);
507 
508  CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
509  vol_ele_slow_rhs, null, null);
510  CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getBoundaryFEName(),
511  natural_bdry_ele_slow_rhs, null, null);
512 
513 
515 }
516 
518  std::string flux_field) {
520  post_proc->addFieldValuesPostProc(mass_field);
521  post_proc->addFieldValuesPostProc(flux_field);
523 }
524 
527  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
528  monitor_ptr, null, null);
530 }
533  // Create solution vector
536  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
537  // Solve problem
538  double ftime = 1;
539  CHKERR TSSetDM(ts, dm);
540  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
541  CHKERR TSSetSolution(ts, X);
542  CHKERR TSSetFromOptions(ts);
543 
544  if (1) {
545  SNES snes;
546  CHKERR TSGetSNES(ts, &snes);
547  KSP ksp;
548  CHKERR SNESGetKSP(snes, &ksp);
549  PC pc;
550  CHKERR KSPGetPC(ksp, &pc);
551  PetscBool is_pcfs = PETSC_FALSE;
552  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
553  // Set up FIELDSPLIT
554  // Only is user set -pc_type fieldsplit
555  if (is_pcfs == PETSC_TRUE) {
556  IS is_mass, is_flux;
557  const MoFEM::Problem *problem_ptr;
558  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
559  CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
560  problem_ptr->getName(), ROW, "MASS1", 0, 1, &is_mass);
561  CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
562  problem_ptr->getName(), ROW, "FLUX1", 0, 1, &is_flux);
563  // CHKERR ISView(is_flux, PETSC_VIEWER_STDOUT_SELF);
564  // CHKERR ISView(is_mass, PETSC_VIEWER_STDOUT_SELF);
565 
566  CHKERR PCFieldSplitSetIS(pc, NULL, is_flux);
567  CHKERR PCFieldSplitSetIS(pc, NULL, is_mass);
568 
569 
570  CHKERR ISDestroy(&is_flux);
571  CHKERR ISDestroy(&is_mass);
572  }
573  }
574 
575  CHKERR TSSolve(ts, X);
577 }
578 
581  global_error = 0;
582  // set nb_species
583  CHKERR setup_system(); // only once
584  nb_species = nb_sp;
585  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
586  CHKERR add_fe("MASS1", "FLUX1"); // nb_species times
587  if (nb_species == 2 || nb_species == 3) {
588  CHKERR add_fe("MASS2", "FLUX2");
589  if (nb_species == 3) {
590  CHKERR add_fe("MASS3", "FLUX3");
591  }
592  }
593  }
594 
595  CHKERR simple_interface->setUp();
596 
597  CHKERR set_blockData(material_blocks);
598 
599  CHKERR extract_bd_ents("ESSENTIAL", "NATURAL", "INTERNAL"); // nb_species times
600 
601  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
602  CHKERR extract_initial_ents(2, inner_surface1);
603  CHKERR update_slow_rhs("MASS1", mass_values_ptr1);
604  if (nb_species == 1) {
605  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
606  "MASS1", data1, data1, data1, material_blocks));
607  } else if (nb_species == 2 || nb_species == 3) {
608  CHKERR extract_initial_ents(3, inner_surface2);
609  CHKERR update_slow_rhs("MASS2", mass_values_ptr2);
610  if (nb_species == 2) {
611  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
612  "MASS1", data1, data2, data2, material_blocks));
613  } else if (nb_species == 3) {
614  CHKERR extract_initial_ents(4, inner_surface3);
615  CHKERR update_slow_rhs("MASS3", mass_values_ptr3);
616  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
617  "MASS1", data1, data2, data3, material_blocks));
618  }
619  }
620  }
621  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
623 
624  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
625  CHKERR push_slow_rhs("MASS1", "FLUX1", data1); // nb_species times
626  if (nb_species == 2 || nb_species == 3) {
627  CHKERR push_slow_rhs("MASS2", "FLUX2", data2);
628  if (nb_species == 3) {
629  CHKERR push_slow_rhs("MASS3", "FLUX3", data3);
630  }
631  }
632  }
633 
634  CHKERR update_vol_fe(vol_ele_stiff_rhs, data1);
635 
636  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
637  CHKERR update_stiff_rhs("MASS1", "FLUX1", mass_values_ptr1,
638  flux_values_ptr1, mass_dots_ptr1, flux_divs_ptr1);
639  CHKERR push_stiff_rhs("MASS1", "FLUX1", data1,
640  material_blocks); // nb_species times
641  if (nb_species == 2 || nb_species == 3) {
642  CHKERR update_stiff_rhs("MASS2", "FLUX2", mass_values_ptr2,
643  flux_values_ptr2, mass_dots_ptr2, flux_divs_ptr2);
644  CHKERR push_stiff_rhs("MASS2", "FLUX2", data2, material_blocks);
645  if (nb_species == 3) {
646  CHKERR update_stiff_rhs("MASS3", "FLUX3", mass_values_ptr3,
647  flux_values_ptr3, mass_dots_ptr3,
648  flux_divs_ptr3);
649  CHKERR push_stiff_rhs("MASS3", "FLUX3", data3, material_blocks);
650  }
651  }
652  }
653 
654  CHKERR update_vol_fe(vol_ele_stiff_lhs, data1);
655 
656  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
657  CHKERR update_stiff_lhs("MASS1", "FLUX1", mass_values_ptr1,
658  flux_values_ptr1);
659  CHKERR push_stiff_lhs("MASS1", "FLUX1", data1,
660  material_blocks); // nb_species times
661  if (nb_species == 2 || nb_species == 3) {
662  CHKERR update_stiff_lhs("MASS2", "FLUX2", mass_values_ptr2,
663  flux_values_ptr2);
664  CHKERR push_stiff_lhs("MASS2", "FLUX2", data2, material_blocks);
665  if (nb_species == 3) {
666  CHKERR update_stiff_lhs("MASS3", "FLUX3", mass_values_ptr3,
667  flux_values_ptr3);
668  CHKERR push_stiff_lhs("MASS3", "FLUX3", data3, material_blocks);
669  }
670  }
671  }
672 
673 
674  // vol_ele_stiff_lhs->getOpPtrVector().push_back(
675  // new OpError(ExactFunction(), ExactFunctionLap(), ExactFunctionGrad(), data1, material_blocks, global_error));
676  CHKERR set_integration_rule();
677  dm = simple_interface->getDM();
678  ts = createTS(m_field.get_comm());
679  boost::shared_ptr<FaceEle> initial_mass_ele(new FaceEle(m_field));
680 
681  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
682  CHKERR apply_IC("MASS1", inner_surface1,
683  initial_mass_ele); // nb_species times
684  if (nb_species == 2 || nb_species == 3) {
685  CHKERR apply_IC("MASS2", inner_surface2, initial_mass_ele);
686  if (nb_species == 3) {
687  CHKERR apply_IC("MASS3", inner_surface3, initial_mass_ele);
688  }
689  }
690  }
691  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
692  initial_mass_ele);
693 
694  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
695  CHKERR apply_BC("FLUX1"); // nb_species times
696  if (nb_species == 2 || nb_species == 3) {
697  CHKERR apply_BC("FLUX2");
698  if (nb_species == 3) {
699  CHKERR apply_BC("FLUX3");
700  }
701  }
702  }
703 
704  CHKERR loop_fe(); // only once
705  post_proc->generateReferenceElementMesh(); // only once
706 
707 
708 
709  if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
710  CHKERR post_proc_fields("MASS1", "FLUX1");
711  post_proc->addFieldValuesPostProc("ERROR");
712  if (nb_species == 2 || nb_species == 3) {
713  CHKERR post_proc_fields("MASS2", "FLUX2");
714  if (nb_species == 3) {
715  CHKERR post_proc_fields("MASS3", "FLUX3");
716  }
717  }
718  }
719 
720  //the double global_vector is a global_vector which sums errors in each element
721 
722  // Vec gVec;
723  // CHKERR VecCreateMPI(PETSC_COMM_WORLD, 1, 1, &gVec);
724  // double err_per_step = 0;
725 
726  // auto compute_global_error = [&gVec](double &g_err, double &err_out) {
727  // MoFEMFunctionBegin;
728 
729  // CHKERR VecZeroEntries(gVec);
730  // CHKERR VecAssemblyBegin(gVec);
731  // CHKERR VecAssemblyEnd(gVec);
732 
733  // int ind[1] = {0};
734  // CHKERR VecSetValues(gVec, 1, ind, &err_out, ADD_VALUES);
735  // CHKERR VecAssemblyBegin(gVec);
736  // CHKERR VecAssemblyEnd(gVec);
737 
738  // CHKERR VecGetValues(gVec, 1, ind, &err_out);
739 
740  // MoFEMFunctionReturn(0);
741  // };
742 
743  // Vec error_per_proc;
744  // CHKERR VecCreateMPI(cOmm, 1, PETSC_DECIDE, &error_per_proc);
745  // auto get_global_error = [&]() {
746  // MoFEMFunctionBegin;
747  // CHKERR VecSetValue(error_per_proc, m_field.get_comm_rank(), global_error, INSERT_VALUES);
748  // MoFEMFunctionReturn(0);
749  // };
750 
751  // auto reinitialize_global = [&]() {
752  // MoFEMFunctionBegin;
753  // CHKERR get_global_error();
754  // CHKERR VecAssemblyBegin(error_per_proc);
755  // CHKERR VecAssemblyEnd(error_per_proc);
756  // double error_sum;
757  // CHKERR VecSum(error_per_proc, &error_sum);
758  // CHKERR PetscPrintf(PETSC_COMM_WORLD, "Error : %3.4e \n",
759  // error_sum);
760  // global_error = 0;
761  // MoFEMFunctionReturn(0);
762  // };
763 
764  // vol_ele_stiff_lhs->postProcessHook = reinitialize_global;
765 
766 
767 
768  monitor_ptr = boost::shared_ptr<Monitor>(
769  new Monitor(cOmm, rAnk, dm, post_proc, global_error)); // nb_species times
770  CHKERR output_result(); // only once
771  CHKERR solve(); // only once
773 }
774 
775 int main(int argc, char *argv[]) {
776  const char param_file[] = "param_file.petsc";
777  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
778  try {
779  moab::Core mb_instance;
780  moab::Interface &moab = mb_instance;
781  MoFEM::Core core(moab);
782  DMType dm_name = "DMMOFEM";
783  CHKERR DMRegister_MoFEM(dm_name);
784 
785  int order = 1;
786  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
787  int nb_species = 1;
788  RDProblem reac_diff_problem(core, order+1);
789  CHKERR reac_diff_problem.run_analysis(nb_species);
790  }
791  CATCH_ERRORS;
793  return 0;
794 }
static Index< 'p', 3 > p
std::string param_file
MoFEM::FaceElementForcesAndSourcesCore FaceEle
@ ROW
Definition: definitions.h:192
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:158
@ L2
field with C-1 continuity
Definition: definitions.h:180
@ HCURL
field with continuous tangents
Definition: definitions.h:178
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ BLOCKSET
Definition: definitions.h:217
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
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
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
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
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:348
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
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 DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1026
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
#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.
int main(int argc, char *argv[])
static char help[]
const double sml
const double ramp_t
const double T
FaceElementForcesAndSourcesCore BoundaryEle
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
OpCalculateJacForFaceImpl< 2 > OpCalculateJacForFace
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
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
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createTS
OpSetContravariantPiolaTransformOnEdge OpSetContrariantPiolaTransformOnEdge
OpCalculateScalarFieldValuesDot OpCalculateScalarValuesDot
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
CoreTmp< 0 > Core
Definition: Core.hpp:1129
OpSetContravariantPiolaTransformFaceImpl< 2 > OpSetContravariantPiolaTransformFace
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
double operator()(const double x) const
double operator()(const double x, const double y, const double t) const
Exact gradient.
FTensor::Tensor1< double, 3 > operator()(const double x, const double y, const double t) const
double operator()(const double x, const double y, const double t) const
double operator()(const double x, const double y, const double t) const
Core (interface) class.
Definition: Core.hpp:77
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
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Deprecated interface functions.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:36
Interface for managing meshsets containing materials and boundary conditions.
Calculate divergence of vector field.
Get vector field for H-div approximation.
Get value at integration points for scalar field.
Make Hdiv space from Hcurl space in 2d.
keeps basic data about problem
std::string getName() const
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
Postprocess on face.
MoFEMErrorCode update_slow_rhs(std::string mass_fiedl, boost::shared_ptr< VectorDouble > &mass_ptr)
boost::shared_ptr< VectorDouble > flux_divs_ptr3
MoFEMErrorCode update_stiff_rhs()
MoFEMErrorCode apply_BC(std::string flux_field)
MoFEMErrorCode add_fe()
MoFEMErrorCode push_slow_rhs()
boost::shared_ptr< VectorDouble > mass_values_ptr3
MoFEMErrorCode update_vol_fe(boost::shared_ptr< FaceEle > &vol_ele, boost::shared_ptr< PreviousData > &data)
MoFEMErrorCode post_proc_fields()
MoFEMErrorCode loop_fe()
MoFEMErrorCode run_analysis()
boost::shared_ptr< MatrixDouble > flux_values_ptr3
MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL)
Range internal_edge_ents
MoFEMErrorCode set_blockData(std::map< int, BlockData > &block_data_map)
RDProblem(MoFEM::Core &core, const int order)
MoFEMErrorCode push_stiff_lhs()
MoFEMErrorCode push_stiff_rhs()
MoFEMErrorCode extract_initial_ents(int block_id, Range &surface)
MoFEMErrorCode solve()
boost::shared_ptr< VectorDouble > mass_dots_ptr3
MoFEMErrorCode apply_IC(std::string mass_field, Range &surface, boost::shared_ptr< FaceEle > &initial_ele, double &init_val)
MoFEMErrorCode update_stiff_lhs(std::string mass_fiedl, std::string flux_field, boost::shared_ptr< VectorDouble > &mass_ptr, boost::shared_ptr< MatrixDouble > &flux_ptr)
MoFEMErrorCode setup_system()
MoFEMErrorCode output_result()
MoFEMErrorCode update_stiff_lhs()
boost::shared_ptr< PreviousData > data3
MoFEMErrorCode update_stiff_rhs(std::string mass_field, std::string flux_field, boost::shared_ptr< VectorDouble > &mass_ptr, boost::shared_ptr< MatrixDouble > &flux_ptr, boost::shared_ptr< VectorDouble > &mass_dot_ptr, boost::shared_ptr< VectorDouble > &flux_div_ptr)
MoFEMErrorCode set_integration_rule()