v0.10.0
elec_phys_2D.cpp
Go to the documentation of this file.
1 #include <stdlib.h>
4 
5 using namespace MoFEM;
6 using namespace ElectroPhysiology;
7 
8 static char help[] = "...\n\n";
9 
10 // #define M_PI 3.14159265358979323846 /* pi */
11 
12 const int dim = 2;
13 
14 double init_val_u = 0.4;
15 double init_val_v = 0.0;
16 
17 // double alpha = -0.08;
18 // double gma = 3;
19 // double ep = 0.005;
20 
21 
22 
23 struct Istimulus {
24  double operator()(const double x, const double y, const double z,
25  const double t) const {
26  if (y <= -1.6000 && t <= 0.50000) {
27  return 80;
28  } else {
29  return 0;
30  }
31  }
32 };
33 
34 struct RhsU {
35  double operator()(const double u, const double v) const {
36  return factor * (1.0 * (c * u * (u - alpha) * (1.0 - u) - u * v));
37  }
38 };
39 
40 struct DRhsU_u {
41  double operator()(const double u, const double v) const {
42  return factor * (c * ((u - alpha) * (1.0 - u) + u * (1.0 - u) - u * (u - alpha)) - v);
43  }
44 };
45 
46 struct RhsV {
47  double operator()(const double u, const double v) const {
48  return factor * ((gma + mu1 * v / (mu2 + u)) * (-v - c * u * (u - b - 1.0)));
49  }
50 };
51 
52 struct DRhsV_v {
53  double operator()(const double u, const double v) const {
54  return factor * (((gma + mu1 / (mu2 + u)) * (-v - c * u * (u - b - 1.0)) -
55  - (gma + mu1 * v / (mu2 + u))));
56  }
57 };
58 
59 const double Tol = 1e-9;
60 
61 struct RK4 {
64  double operator()(const double u, const double v, const double dt) const {
65 
66 
67  double dv;
68 
69  double vk = v;
70 
71  double err = 1;
72 
73  while (err > Tol){
74  double rhsv = rhs_v(u, vk);
75  double Drhsv = Drhs_v(u, vk);
76  dv = - 1.0 / (1 - dt * Drhsv) * (vk - v - dt * rhsv);
77  vk += dv;
78  err = abs(dv);
79  }
80 
81  return vk;
82 
83  // double k1 = dt * rhs_v(u, v);
84  // // double k2 = dt * rhs_v(u, v + 0.5 * k1);
85  // // double k3 = dt * rhs_v(u, v + 0.5 * k2);
86  // double k4 = dt * rhs_v(u, v + k1);
87  // return v + 0.5 * (k1 + k4);
88  }
89 };
90 
91 struct RDProblem {
92 public:
93  RDProblem(MoFEM::Core &core, const int order)
94  : m_field(core), order(order), cOmm(m_field.get_comm()),
95  rAnk(m_field.get_comm_rank()) {
96  vol_ele_slow_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
97  natural_bdry_ele_slow_rhs =
98  boost::shared_ptr<BoundaryEle>(new BoundaryEle(m_field));
99  vol_ele_stiff_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
100  vol_ele_stiff_lhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
101  post_proc = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
102  new PostProcFaceOnRefinedMesh(m_field));
103 
104  data1 = boost::shared_ptr<PreviousData>(new PreviousData());
105  data2 = boost::shared_ptr<PreviousData>(new PreviousData());
106 
107 
108  flux_values_ptr1 =
109  boost::shared_ptr<MatrixDouble>(data1, &data1->flux_values);
110 
111  flux_divs_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->flux_divs);
112 
113  mass_values_ptr1 =
114  boost::shared_ptr<VectorDouble>(data1, &data1->mass_values);
115 
116  mass_dots_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->mass_dots);
117 
118  flux_values_ptr2 =
119  boost::shared_ptr<MatrixDouble>(data2, &data2->flux_values);
120  flux_divs_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->flux_divs);
121 
122  mass_values_ptr2 =
123  boost::shared_ptr<VectorDouble>(data2, &data2->mass_values);
124  mass_dots_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->mass_dots);
125 
126 
127  }
128 
129  // RDProblem(const int order) : order(order){}
130  MoFEMErrorCode run_analysis();
131 
132  double global_error;
133 
134 private:
135  MoFEMErrorCode setup_system();
136  MoFEMErrorCode add_fe();
137  MoFEMErrorCode set_blockData(std::map<int, BlockData> &block_data_map);
138  MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL);
139  MoFEMErrorCode extract_initial_ents(int block_id, Range &surface);
140  MoFEMErrorCode update_slow_rhs(std::string mass_fiedl,
141  boost::shared_ptr<VectorDouble> &mass_ptr);
142  MoFEMErrorCode push_slow_rhs(std::string mass_field, std::string flux_field,
143  boost::shared_ptr<PreviousData> &data1,
144  boost::shared_ptr<PreviousData> &data2);
145  MoFEMErrorCode update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
146  boost::shared_ptr<PreviousData> &data);
148  update_stiff_rhs(std::string mass_field, std::string flux_field,
149  boost::shared_ptr<VectorDouble> &mass_ptr,
150  boost::shared_ptr<MatrixDouble> &flux_ptr,
151  boost::shared_ptr<VectorDouble> &mass_dot_ptr,
152  boost::shared_ptr<VectorDouble> &flux_div_ptr);
153  MoFEMErrorCode push_stiff_rhs(std::string mass_field, std::string flux_field,
154  boost::shared_ptr<PreviousData> &data1,
155  boost::shared_ptr<PreviousData> &data2,
156  std::map<int, BlockData> &block_map,
157  Range &surface);
158  MoFEMErrorCode update_stiff_lhs(std::string mass_fiedl,
159  std::string flux_field,
160  boost::shared_ptr<VectorDouble> &mass_ptr,
161  boost::shared_ptr<MatrixDouble> &flux_ptr);
162  MoFEMErrorCode push_stiff_lhs(std::string mass_field, std::string flux_field,
163  boost::shared_ptr<PreviousData> &datau,
164  boost::shared_ptr<PreviousData> &datav,
165  std::map<int, BlockData> &block_map);
166 
167  MoFEMErrorCode set_integration_rule();
168  MoFEMErrorCode apply_IC(std::string mass_field, Range &surface,
169  boost::shared_ptr<FaceEle> &initial_ele,
170  double & init_val);
171  MoFEMErrorCode apply_BC(std::string flux_field);
172  MoFEMErrorCode loop_fe();
173  MoFEMErrorCode post_proc_fields();
174  MoFEMErrorCode output_result();
175  MoFEMErrorCode solve();
176 
181 
184 
185  Range inner_surface1; // nb_species times
187 
188 
189  MPI_Comm cOmm;
190  const int rAnk;
191 
192  int order;
194 
195  std::map<int, BlockData> material_blocks;
196 
197  boost::shared_ptr<FaceEle> vol_ele_slow_rhs;
198  boost::shared_ptr<FaceEle> vol_ele_stiff_rhs;
199  boost::shared_ptr<FaceEle> vol_ele_stiff_lhs;
200  boost::shared_ptr<BoundaryEle> natural_bdry_ele_slow_rhs;
201 
202  boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc;
203  boost::shared_ptr<Monitor> monitor_ptr;
204 
205  boost::shared_ptr<PreviousData> data1; // nb_species times
206  boost::shared_ptr<PreviousData> data2;
207 
208 
209  boost::shared_ptr<MatrixDouble> flux_values_ptr1; // nb_species times
210  boost::shared_ptr<MatrixDouble> flux_values_ptr2;
211 
212 
213  boost::shared_ptr<VectorDouble> flux_divs_ptr1; // nb_species times
214  boost::shared_ptr<VectorDouble> flux_divs_ptr2;
215 
216 
217  boost::shared_ptr<VectorDouble> mass_values_ptr1; // nb_species times
218  boost::shared_ptr<VectorDouble> mass_values_ptr2;
219 
220 
221  boost::shared_ptr<VectorDouble> mass_dots_ptr1; // nb_species times
222  boost::shared_ptr<VectorDouble> mass_dots_ptr2;
223 
224 
225  boost::shared_ptr<ForcesAndSourcesCore> null;
226 };
227 
230  CHKERR m_field.getInterface(simple_interface);
231  CHKERR simple_interface->getOptions();
232  CHKERR simple_interface->loadFile();
234 }
235 
238  CHKERR simple_interface->addDomainField("U", L2,
240  CHKERR simple_interface->addDataField("V", L2,
242 
243  CHKERR simple_interface->addDomainField("F", HCURL,
245 
246  CHKERR simple_interface->addBoundaryField("F", HCURL,
248 
249 
250 
251  CHKERR simple_interface->setFieldOrder("U", order - 1);
252  CHKERR simple_interface->setFieldOrder("V", order - 1);
253  CHKERR simple_interface->setFieldOrder("F", order);
254 
255 
257 }
258 
262  string name = it->getName();
263  const int id = it->getMeshsetId();
264  if (name.compare(0, 14, "REGION1") == 0) {
265  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
266  id, BLOCKSET, 2, block_map[id].block_ents, true);
267  block_map[id].B0 = 1e-3;
268  block_map[id].block_id = id;
269  } else if (name.compare(0, 14, "REGION2") == 0) {
270  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
271  id, BLOCKSET, 2, block_map[id].block_ents, true);
272  block_map[id].B0 = 1e-1;
273  block_map[id].block_id = id;
274  }
275  }
277 }
278 
280  std::string natural) {
283  string name = it->getName();
284  if (name.compare(0, 14, natural) == 0) {
285 
286  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
287  natural_bdry_ents, true);
288  } else if (name.compare(0, 14, essential) == 0) {
289  CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
290  essential_bdry_ents, true);
291  }
292  }
294 }
295 
296 MoFEMErrorCode RDProblem::extract_initial_ents(int block_id, Range &surface) {
298  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(block_id,
299  BLOCKSET)) {
300  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
301  block_id, BLOCKSET, 2, surface, true);
302  }
304 }
305 
307 RDProblem::update_slow_rhs(std::string mass_field,
308  boost::shared_ptr<VectorDouble> &mass_ptr) {
310  vol_ele_slow_rhs->getOpPtrVector().push_back(
311  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
313 }
314 
316 RDProblem::push_slow_rhs(std::string mass_field, std::string flux_field,
317  boost::shared_ptr<PreviousData> &data_1,
318  boost::shared_ptr<PreviousData> &data_2) {
320 
321  vol_ele_slow_rhs->getOpPtrVector().push_back(
322  new OpAssembleSlowRhsV(mass_field, data_1, data_2, RhsU()));
323 
324  // natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
325  // new OpAssembleNaturalBCRhsTau(flux_field, natural_bdry_ents));
326 
327  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
328  new OpEssentialBC(flux_field, essential_bdry_ents));
329 
331 }
332 
333 MoFEMErrorCode RDProblem::update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
334  boost::shared_ptr<PreviousData> &data) {
336  vol_ele->getOpPtrVector().push_back(new OpCalculateJacForFace(data->jac));
337  vol_ele->getOpPtrVector().push_back(
338  new OpCalculateInvJacForFace(data->inv_jac));
339  vol_ele->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
340 
341  vol_ele->getOpPtrVector().push_back(
342  new OpSetContravariantPiolaTransformFace(data->jac));
343 
344  vol_ele->getOpPtrVector().push_back(new OpSetInvJacHcurlFace(data->inv_jac));
346 }
347 
349 RDProblem::update_stiff_rhs(std::string mass_field, std::string flux_field,
350  boost::shared_ptr<VectorDouble> &mass_ptr,
351  boost::shared_ptr<MatrixDouble> &flux_ptr,
352  boost::shared_ptr<VectorDouble> &mass_dot_ptr,
353  boost::shared_ptr<VectorDouble> &flux_div_ptr) {
354 
356 
357  vol_ele_stiff_rhs->getOpPtrVector().push_back(
358  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
359 
360  vol_ele_stiff_rhs->getOpPtrVector().push_back(
361  new OpCalculateHdivVectorField<3>(flux_field, flux_ptr));
362 
363  vol_ele_stiff_rhs->getOpPtrVector().push_back(
364  new OpCalculateScalarValuesDot(mass_field, mass_dot_ptr));
365 
366  vol_ele_stiff_rhs->getOpPtrVector().push_back(
367  new OpCalculateHdivVectorDivergence<3, 2>(flux_field, flux_div_ptr));
369 }
370 
372  std::string flux_field,
373  boost::shared_ptr<PreviousData> &data1,
374  boost::shared_ptr<PreviousData> &data2,
375  std::map<int, BlockData> &block_map,
376  Range &surface) {
378  vol_ele_stiff_rhs->getOpPtrVector().push_back(
379  new OpAssembleStiffRhsTau<3>(flux_field, data1, block_map));
380 
381  vol_ele_stiff_rhs->getOpPtrVector().push_back(
382  new OpAssembleStiffRhsV<3>(mass_field, data1, data2, RhsU(), block_map, surface));
384 }
385 
387 RDProblem::update_stiff_lhs(std::string mass_field, std::string flux_field,
388  boost::shared_ptr<VectorDouble> &mass_ptr,
389  boost::shared_ptr<MatrixDouble> &flux_ptr) {
391  vol_ele_stiff_lhs->getOpPtrVector().push_back(
392  new OpCalculateScalarFieldValues(mass_field, mass_ptr));
393 
394  vol_ele_stiff_lhs->getOpPtrVector().push_back(
395  new OpCalculateHdivVectorField<3>(flux_field, flux_ptr));
397 }
398 
400  std::string flux_field,
401  boost::shared_ptr<PreviousData> &data1,
402  boost::shared_ptr<PreviousData> &data2,
403  std::map<int, BlockData> &block_map) {
405  vol_ele_stiff_lhs->getOpPtrVector().push_back(
406  new OpAssembleLhsTauTau<3>(flux_field, data1, block_map));
407 
408  vol_ele_stiff_lhs->getOpPtrVector().push_back(
409  new OpAssembleLhsVV(mass_field, data1, data2, DRhsU_u()));
410 
411  vol_ele_stiff_lhs->getOpPtrVector().push_back(
412  new OpAssembleLhsTauV<3>(flux_field, mass_field, data1, block_map));
413 
414  vol_ele_stiff_lhs->getOpPtrVector().push_back(
415  new OpAssembleLhsVTau(mass_field, flux_field));
417 }
418 
419 
422  auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
423  vol_ele_slow_rhs->getRuleHook = vol_rule;
424  natural_bdry_ele_slow_rhs->getRuleHook = vol_rule;
425 
426  vol_ele_stiff_rhs->getRuleHook = vol_rule;
427 
428  vol_ele_stiff_lhs->getRuleHook = vol_rule;
430 }
431 
432 MoFEMErrorCode RDProblem::apply_IC(std::string mass_field, Range &surface,
433  boost::shared_ptr<FaceEle> &initial_ele,
434  double & init_val) {
436  initial_ele->getOpPtrVector().push_back(
437  new OpInitialMass(mass_field, surface, init_val));
439 }
440 
441 MoFEMErrorCode RDProblem::apply_BC(std::string flux_field) {
443  CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
444  "SimpleProblem", flux_field, essential_bdry_ents);
445 
447 }
450  CHKERR TSSetType(ts, TSARKIMEX);
451  CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
452 
453  CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
454  vol_ele_stiff_lhs, null, null);
455 
456  CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
457  vol_ele_stiff_rhs, null, null);
458 
459  CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
460  vol_ele_slow_rhs, null, null);
461  CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getBoundaryFEName(),
462  natural_bdry_ele_slow_rhs, null, null);
464 }
465 
468  post_proc->addFieldValuesPostProc("U");
469  post_proc->addFieldValuesPostProc("F");
470  post_proc->addFieldValuesPostProc("V");
471 
473  }
474 
477  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
478  monitor_ptr, null, null);
480  }
483  // Create solution vector
486  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
487  // Solve problem
488  double ftime = 1;
489  CHKERR TSSetDM(ts, dm);
490  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
491  CHKERR TSSetSolution(ts, X);
492  CHKERR TSSetFromOptions(ts);
493 
494  if (1) {
495  SNES snes;
496  CHKERR TSGetSNES(ts, &snes);
497  KSP ksp;
498  CHKERR SNESGetKSP(snes, &ksp);
499  PC pc;
500  CHKERR KSPGetPC(ksp, &pc);
501  PetscBool is_pcfs = PETSC_FALSE;
502  PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
503  // Set up FIELDSPLIT
504  // Only is user set -pc_type fieldsplit
505  if (is_pcfs == PETSC_TRUE) {
506  IS is_mass, is_flux;
507  const MoFEM::Problem *problem_ptr;
508  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
509  CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
510  problem_ptr->getName(), ROW, "U", 0, 1, &is_mass);
511  CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
512  problem_ptr->getName(), ROW, "F", 0, 1, &is_flux);
513  // CHKERR ISView(is_flux, PETSC_VIEWER_STDOUT_SELF);
514  // CHKERR ISView(is_mass, PETSC_VIEWER_STDOUT_SELF);
515 
516  CHKERR PCFieldSplitSetIS(pc, NULL, is_mass);
517  CHKERR PCFieldSplitSetIS(pc, NULL, is_flux);
518 
519 
520  CHKERR ISDestroy(&is_flux);
521  CHKERR ISDestroy(&is_mass);
522  }
523  }
524 
525  CHKERR TSSolve(ts, X);
527  }
528 
531  global_error = 0;
532 
533  CHKERR setup_system(); // only once
534 
535  CHKERR add_fe(); // nb_species times
536 
537  CHKERR simple_interface->setUp();
538 
539  CHKERR set_blockData(material_blocks);
540 
541  CHKERR extract_bd_ents("ESSENTIAL", "NATURAL"); // nb_species times
542 
543  CHKERR extract_initial_ents(2, inner_surface1);
544  CHKERR extract_initial_ents(3, inner_surface2);
545 
546  CHKERR update_slow_rhs("U", mass_values_ptr1);
547  CHKERR update_slow_rhs("V", mass_values_ptr2);
548 
549  natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
551 
552  CHKERR push_slow_rhs("U", "F", data1, data2); // nb_species times
553 
554  CHKERR update_vol_fe(vol_ele_stiff_rhs, data1);
555 
556  CHKERR update_stiff_rhs("U", "F", mass_values_ptr1,
557  flux_values_ptr1, mass_dots_ptr1, flux_divs_ptr1);
558  vol_ele_stiff_rhs->getOpPtrVector().push_back(
559  new OpCalculateScalarFieldValues("V", mass_values_ptr2));
560  CHKERR push_stiff_rhs("U", "F", data1, data2,
561  material_blocks, inner_surface2); // nb_species times
562 
563  CHKERR update_vol_fe(vol_ele_stiff_lhs, data1);
564 
565  CHKERR update_stiff_lhs("U", "F", mass_values_ptr1,
566  flux_values_ptr1);
567 
568  vol_ele_stiff_lhs->getOpPtrVector().push_back(
569  new OpCalculateScalarFieldValues("V", mass_values_ptr2));
570  CHKERR push_stiff_lhs("U", "F", data1, data2,
571  material_blocks); // nb_species times
572  vol_ele_stiff_lhs->getOpPtrVector().push_back(
573  new OpSolveRecovery("V", data1, data2, RK4()));
574 
575  CHKERR set_integration_rule();
576  dm = simple_interface->getDM();
577  ts = createTS(m_field.get_comm());
578  boost::shared_ptr<FaceEle> initial_mass_ele(new FaceEle(m_field));
579 
580  CHKERR apply_IC("U", inner_surface1, initial_mass_ele,
581  init_val_u); // nb_species times
582  CHKERR apply_IC("V", inner_surface1, initial_mass_ele,
583  init_val_v); // nb_species times
584 
585  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
586  initial_mass_ele);
587 
588  CHKERR apply_BC("F"); // nb_species times
589 
590  CHKERR loop_fe(); // only once
591  post_proc->generateReferenceElementMesh(); // only once
592 
593  CHKERR post_proc_fields();
594 
595  monitor_ptr = boost::shared_ptr<Monitor>(
596  new Monitor(cOmm, rAnk, dm, post_proc)); // nb_species times
597  CHKERR output_result(); // only once
598  CHKERR solve(); // only once
600  }
601 
602  int main(int argc, char *argv[]) {
603  const char param_file[] = "param_file.petsc";
604  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
605  try {
606  moab::Core mb_instance;
607  moab::Interface &moab = mb_instance;
608  MoFEM::Core core(moab);
609  DMType dm_name = "DMMOFEM";
610  CHKERR DMRegister_MoFEM(dm_name);
611 
612  int order = 1;
613  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
614  RDProblem reac_diff_problem(core, order + 1);
615  CHKERR reac_diff_problem.run_analysis();
616  }
617  CATCH_ERRORS;
619  return 0;
620  }
boost::shared_ptr< FaceEle > vol_ele_stiff_lhs
boost::shared_ptr< VectorDouble > flux_divs_ptr2
Calculate divergence of vector field.
Range natural_bdry_ents
Range essential_bdry_ents
boost::shared_ptr< MatrixDouble > flux_values_ptr1
Deprecated interface functions.
Problem manager is used to build and partition problems.
MoFEMErrorCode add_fe()
double operator()(const double u, const double v) const
MoFEMErrorCode update_stiff_lhs()
CoreTmp< 0 > Core
Definition: Core.hpp:1129
MoFEMErrorCode setup_system()
const double Tol
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
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
MPI_Comm cOmm
boost::shared_ptr< VectorDouble > flux_divs_ptr1
double init_val_u
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 char help[]
Definition: elec_phys_2D.cpp:8
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
Range inner_surface2
boost::shared_ptr< PostProcFaceOnRefinedMesh > post_proc
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,...
boost::shared_ptr< Monitor > monitor_ptr
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
boost::shared_ptr< FaceEle > vol_ele_stiff_rhs
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
MoFEMErrorCode apply_BC(std::string flux_field)
SmartPetscObj< DM > dm
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
DRhsV_v Drhs_v
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)
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
const int dim
MoFEM::Interface & m_field
RhsV rhs_v
Range inner_surface1
boost::shared_ptr< VectorDouble > mass_values_ptr2
boost::shared_ptr< FaceEle > vol_ele_slow_rhs
MoFEMErrorCode push_stiff_rhs()
RDProblem(MoFEM::Core &core, const int order)
constexpr int order
MoFEMErrorCode extract_initial_ents(int block_id, Range &surface)
Calculate inverse of jacobian for face element.
boost::shared_ptr< VectorDouble > mass_values_ptr1
MoFEM::FaceElementForcesAndSourcesCore FaceEle
double operator()(const double u, const double v) const
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()
Simple * simple_interface
boost::shared_ptr< BoundaryEle > natural_bdry_ele_slow_rhs
std::map< int, BlockData > material_blocks
MoFEMErrorCode loop_fe()
field with continuous tangents
Definition: definitions.h:178
Get vector field for H-div approximation.
const int rAnk
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:604
MoFEMErrorCode update_stiff_rhs()
boost::shared_ptr< MatrixDouble > flux_values_ptr2
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
SmartPetscObj< TS > ts
int main(int argc, char *argv[])
MoFEMErrorCode output_result()
boost::shared_ptr< VectorDouble > mass_dots_ptr2
Postprocess on face.
MoFEMErrorCode solve()
brief Transform local reference derivatives of shape function to global derivatives for face
double operator()(const double u, const double v) const
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
double init_val_v
boost::shared_ptr< PreviousData > data2
std::string param_file
double operator()(const double u, const double v, const double dt) const
Get value at integration points for scalar field.
OpCalculateScalarFieldValuesDot OpCalculateScalarValuesDot
boost::shared_ptr< VectorDouble > mass_dots_ptr1
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
double global_error
double operator()(const double u, const double v) const
double operator()(const double x, const double y, const double z, const double t) const
Make Hdiv space from Hcurl space in 2d.
boost::shared_ptr< PreviousData > data1
field with C-1 continuity
Definition: definitions.h:180