v0.14.0
reaction_diffusion.cpp
Go to the documentation of this file.
1 /**
2  * \file reaction_diffusion.cpp
3  * \example reaction_diffusion.cpp
4  *
5  **/
6 
8 using namespace MoFEM;
9 
10 static char help[] = "...\n\n";
11 
13 
17 
19 
20 const double D = 2e-3; ///< diffusivity
21 const double r = 1; ///< rate factor
22 const double k = 1; ///< caring capacity
23 
24 const double u0 = 0.1; ///< inital vale on blocksets
25 
26 const int save_every_nth_step = 1;
27 
28 /**
29  * @brief Common data
30  *
31  * Common data are used to keep and pass data between elements
32  *
33  */
34 struct CommonData {
35 
36  MatrixDouble grad; ///< Gradients of field "u" at integration points
37  VectorDouble val; ///< Values of field "u" at integration points
38  VectorDouble dot_val; ///< Rate of values of field "u" at integration points
39 
40  SmartPetscObj<Mat> M; ///< Mass matrix
41  SmartPetscObj<KSP> ksp; ///< Linear solver
42 };
43 
44 /**
45  * @brief Assemble mass matrix
46  */
48  OpAssembleMass(boost::shared_ptr<CommonData> &data)
49  : OpEle("u", "u", OpEle::OPROWCOL), commonData(data) {
50  sYmm = true;
51  }
52  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
53  EntityType col_type, EntData &row_data,
54  EntData &col_data) {
56  const int nb_row_dofs = row_data.getIndices().size();
57  const int nb_col_dofs = col_data.getIndices().size();
58  if (nb_row_dofs && nb_col_dofs) {
59  const int nb_integration_pts = getGaussPts().size2();
60  mat.resize(nb_row_dofs, nb_col_dofs, false);
61  mat.clear();
62  auto t_row_base = row_data.getFTensor0N();
63  auto t_w = getFTensor0IntegrationWeight();
64  const double vol = getMeasure();
65  for (int gg = 0; gg != nb_integration_pts; ++gg) {
66  const double a = t_w * vol;
67  for (int rr = 0; rr != nb_row_dofs; ++rr) {
68  auto t_col_base = col_data.getFTensor0N(gg, 0);
69  for (int cc = 0; cc != nb_col_dofs; ++cc) {
70  mat(rr, cc) += a * t_row_base * t_col_base;
71  ++t_col_base;
72  }
73  ++t_row_base;
74  }
75  ++t_w;
76  }
77 
78  CHKERR MatSetValues(commonData->M, row_data, col_data, &mat(0, 0),
79  ADD_VALUES);
80  if (row_side != col_side || row_type != col_type) {
81  transMat.resize(nb_col_dofs, nb_row_dofs, false);
82  noalias(transMat) = trans(mat);
83  CHKERR MatSetValues(commonData->M, col_data, row_data, &transMat(0, 0),
84  ADD_VALUES);
85  }
86  }
88  }
89 
90 private:
92  boost::shared_ptr<CommonData> commonData;
93 };
94 
95 /**
96  * @brief Assemble slow part
97  *
98  * Solve problem \f$ F(t,u,\dot{u}) = G(t,u) \f$ where here the right hand side
99  * \f$ G(t,u) \f$ is implemented.
100  *
101  */
103  OpAssembleSlowRhs(boost::shared_ptr<CommonData> &data)
104  : OpEle("u", OpEle::OPROW), commonData(data) {}
105  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
107  const int nb_dofs = data.getIndices().size();
108  if (nb_dofs) {
109  vecF.resize(nb_dofs, false);
110  vecF.clear();
111 
112  const int nb_integration_pts = getGaussPts().size2();
113  auto t_val = getFTensor0FromVec(commonData->val);
114  auto t_base = data.getFTensor0N();
115  auto t_w = getFTensor0IntegrationWeight();
116 
117  const double vol = getMeasure();
118  for (int gg = 0; gg != nb_integration_pts; ++gg) {
119  const double a = vol * t_w;
120  const double f = a * r * t_val * (1 - t_val / k);
121  for (int rr = 0; rr != nb_dofs; ++rr) {
122  const double b = f * t_base;
123  vecF[rr] += b;
124  ++t_base;
125  }
126 
127  ++t_val;
128  ++t_w;
129  }
130 
131  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
132  PETSC_TRUE);
133  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
134  ADD_VALUES);
135  }
137  }
138 
139 private:
140  boost::shared_ptr<CommonData> commonData;
142 };
143 
144 /**
145  * @brief Assemble stiff part
146  *
147  * Solve problem \f$ F(t,u,\dot{u}) = G(t,u) \f$ where here the right hand side
148  * \f$ F(t,u,\dot{u}) \f$ is implemented.
149  *
150  */
151 template <int DIM> struct OpAssembleStiffRhs : OpEle {
152  OpAssembleStiffRhs(boost::shared_ptr<CommonData> &data)
153  : OpEle("u", OpEle::OPROW), commonData(data) {}
154  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
156  const int nb_dofs = data.getIndices().size();
157  if (nb_dofs) {
158  vecF.resize(nb_dofs, false);
159  vecF.clear();
160 
161  const int nb_integration_pts = getGaussPts().size2();
162  auto t_dot_val = getFTensor0FromVec(commonData->dot_val);
163  auto t_grad = getFTensor1FromMat<DIM>(commonData->grad);
164  auto t_base = data.getFTensor0N();
165  auto t_diff_base = data.getFTensor1DiffN<DIM>();
166  auto t_w = getFTensor0IntegrationWeight();
168 
169  const double vol = getMeasure();
170  for (int gg = 0; gg != nb_integration_pts; ++gg) {
171  const double a = vol * t_w;
172  for (int rr = 0; rr != nb_dofs; ++rr) {
173  vecF[rr] += a * (t_base * t_dot_val + D * t_diff_base(i) * t_grad(i));
174  ++t_diff_base;
175  ++t_base;
176  }
177  ++t_dot_val;
178  ++t_grad;
179  ++t_w;
180  }
181 
182  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
183  PETSC_TRUE);
184  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
185  ADD_VALUES);
186  }
188  }
189 
190 private:
191  boost::shared_ptr<CommonData> commonData;
193 };
194 
195 /**
196  * @brief Assemble stiff part tangent
197  *
198  * Solve problem \f$ F(t,u,\dot{u}) = G(t,u) \f$ where here the right hand side
199  * \f$ \frac{\textrm{d} F}{\textrm{d} u^n} = a F_{\dot{u}}(t,u,\textrm{u}) +
200  * F_{u}(t,u,\textrm{u}) \f$ is implemented.
201  *
202  */
203 template <int DIM> struct OpAssembleStiffLhs : OpEle {
204  OpAssembleStiffLhs(boost::shared_ptr<CommonData> &data)
205  : OpEle("u", "u", OpEle::OPROWCOL), commonData(data) {
206  sYmm = true;
207  }
208  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
209  EntityType col_type, EntData &row_data,
210  EntData &col_data) {
212  const int nb_row_dofs = row_data.getIndices().size();
213  const int nb_col_dofs = col_data.getIndices().size();
214  if (nb_row_dofs && nb_col_dofs) {
215  mat.resize(nb_row_dofs, nb_col_dofs, false);
216  mat.clear();
217 
218  const int nb_integration_pts = getGaussPts().size2();
219  auto t_row_base = row_data.getFTensor0N();
220  auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
221  auto t_w = getFTensor0IntegrationWeight();
223 
224  const double ts_a = getFEMethod()->ts_a;
225  const double vol = getMeasure();
226  for (int gg = 0; gg != nb_integration_pts; ++gg) {
227  const double a = vol * t_w;
228  for (int rr = 0; rr != nb_row_dofs; ++rr) {
229 
230  auto t_col_base = col_data.getFTensor0N(gg, 0);
231  auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
232 
233  for (int cc = 0; cc != nb_col_dofs; ++cc) {
234  mat(rr, cc) += a * (t_row_base * t_col_base * ts_a +
235  D * t_row_diff_base(i) * t_col_diff_base(i));
236  ++t_col_base;
237  ++t_col_diff_base;
238  }
239 
240  ++t_row_base;
241  ++t_row_diff_base;
242  }
243  ++t_w;
244  }
245 
246  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
247  ADD_VALUES);
248  if (row_side != col_side || row_type != col_type) {
249  transMat.resize(nb_col_dofs, nb_row_dofs, false);
250  noalias(transMat) = trans(mat);
251  CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
252  &transMat(0, 0), ADD_VALUES);
253  }
254  }
256  }
257 
258 private:
259  boost::shared_ptr<CommonData> commonData;
261 };
262 
263 /**
264  * @brief Monitor solution
265  *
266  * This functions is called by TS solver at the end of each step. It is used
267  * to output results to the hard drive.
268  */
269 struct Monitor : public FEMethod {
270 
271  Monitor(SmartPetscObj<DM> &dm, boost::shared_ptr<PostProcEle> &post_proc)
272  : dM(dm), postProc(post_proc){};
273 
274  MoFEMErrorCode preProcess() { return 0; }
275  MoFEMErrorCode operator()() { return 0; }
276 
279  if (ts_step % save_every_nth_step == 0) {
280  CHKERR DMoFEMLoopFiniteElements(dM, "dFE", postProc);
281  CHKERR postProc->writeFile(
282  "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
283  }
285  }
286 
287 private:
289  boost::shared_ptr<PostProcEle> postProc;
290 };
291 
292 }; // namespace ReactionDiffusionEquation
293 
294 using namespace ReactionDiffusionEquation;
295 
296 int main(int argc, char *argv[]) {
297 
298  // initialize petsc
299  const char param_file[] = "param_file.petsc";
300  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
301 
302  try {
303 
304  // Create moab and mofem instances
305  moab::Core mb_instance;
306  moab::Interface &moab = mb_instance;
307  MoFEM::Core core(moab);
308  MoFEM::Interface &m_field = core;
309 
310  // Register DM Manager
311  DMType dm_name = "DMMOFEM";
312  CHKERR DMRegister_MoFEM(dm_name);
313 
314  // Simple interface
315  Simple *simple_interface;
316  CHKERR m_field.getInterface(simple_interface);
317  CHKERR simple_interface->getOptions();
318  CHKERR simple_interface->loadFile();
319 
320  int order = 4; ///< approximation order
321  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
322 
323  // add fields
324  CHKERR simple_interface->addDomainField("u", H1, AINSWORTH_LEGENDRE_BASE,
325  1);
326  // set fields order
327  CHKERR simple_interface->setFieldOrder("u", order);
328  // setup problem
329  CHKERR simple_interface->setUp();
330 
331  // Create common data structure
332  boost::shared_ptr<CommonData> data(new CommonData());
333  /// Alias pointers to data in common data structure
334  auto val_ptr = boost::shared_ptr<VectorDouble>(data, &data->val);
335  auto dot_val_ptr = boost::shared_ptr<VectorDouble>(data, &data->dot_val);
336  auto grad_ptr = boost::shared_ptr<MatrixDouble>(data, &data->grad);
337 
338  // Create finite element instances to integrate the right-hand side of slow
339  // and stiff vector, and the tangent left-hand side for stiff part.
340  boost::shared_ptr<Ele> vol_ele_slow_rhs(new Ele(m_field));
341  boost::shared_ptr<Ele> vol_ele_stiff_rhs(new Ele(m_field));
342  boost::shared_ptr<Ele> vol_ele_stiff_lhs(new Ele(m_field));
343 
344  // Push operators to integrate the slow right-hand side vector
345  vol_ele_slow_rhs->getOpPtrVector().push_back(
346  new OpCalculateScalarFieldValues("u", val_ptr));
347  vol_ele_slow_rhs->getOpPtrVector().push_back(new OpAssembleSlowRhs(data));
348 
349  // PETSc IMAX and Explicit solver demans that g = M^-1 G is provided. So
350  // when the slow right-hand side vector (G) is assembled is solved for g
351  // vector.
352  auto solve_for_g = [&]() {
354  if (*(vol_ele_slow_rhs->vecAssembleSwitch)) {
355  CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
356  SCATTER_REVERSE);
357  CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
358  SCATTER_REVERSE);
359  CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
360  CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
361  *vol_ele_slow_rhs->vecAssembleSwitch = false;
362  }
363  CHKERR KSPSolve(data->ksp, vol_ele_slow_rhs->ts_F,
364  vol_ele_slow_rhs->ts_F);
366  };
367  // Add hook to the element to calculate g.
368  vol_ele_slow_rhs->postProcessHook = solve_for_g;
369 
370  auto det_ptr = boost::make_shared<VectorDouble>();
371  auto jac_ptr = boost::make_shared<MatrixDouble>();
372  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
373  // Add operators to calculate the stiff right-hand side
374  vol_ele_stiff_rhs->getOpPtrVector().push_back(
375  new OpCalculateHOJac<2>(jac_ptr));
376  vol_ele_stiff_rhs->getOpPtrVector().push_back(
377  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
378  vol_ele_stiff_rhs->getOpPtrVector().push_back(
379  new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
380  vol_ele_stiff_rhs->getOpPtrVector().push_back(
381  new OpCalculateScalarFieldValuesDot("u", dot_val_ptr));
382  vol_ele_stiff_rhs->getOpPtrVector().push_back(
383  new OpCalculateScalarFieldGradient<2>("u", grad_ptr));
384  vol_ele_stiff_rhs->getOpPtrVector().push_back(
385  new OpAssembleStiffRhs<2>(data));
386 
387  // Add operators to calculate the stiff left-hand side
388  vol_ele_stiff_lhs->getOpPtrVector().push_back(
389  new OpCalculateHOJac<2>(jac_ptr));
390  vol_ele_stiff_lhs->getOpPtrVector().push_back(
391  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
392  vol_ele_stiff_lhs->getOpPtrVector().push_back(
393  new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
394  vol_ele_stiff_lhs->getOpPtrVector().push_back(
395  new OpAssembleStiffLhs<2>(data));
396 
397  // Set integration rules
398  auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
399  vol_ele_slow_rhs->getRuleHook = vol_rule;
400  vol_ele_stiff_rhs->getRuleHook = vol_rule;
401  vol_ele_stiff_lhs->getRuleHook = vol_rule;
402 
403  // Crate element for post-processing
404 
405  auto post_proc = boost::make_shared<PostProcEle>(m_field);
406  boost::shared_ptr<ForcesAndSourcesCore> null;
407 
409 
410  auto u_ptr = boost::make_shared<VectorDouble>();
411  post_proc->getOpPtrVector().push_back(
412  new OpCalculateScalarFieldValues("u", u_ptr));
413  post_proc->getOpPtrVector().push_back(
414 
415  new OpPPMap(
416 
417  post_proc->getPostProcMesh(), post_proc->getMapGaussPts(),
418 
419  {{"u", u_ptr}},
420 
421  {},
422 
423  {},
424 
425  {}
426 
427  )
428 
429  );
430 
431  // Get PETSc discrete manager
432  auto dm = simple_interface->getDM();
433 
434  // Get surface entities form blockset, set initial values in those
435  // blocksets. To keep it simple is assumed that inital values are on
436  // blockset 1
437  if (m_field.getInterface<MeshsetsManager>()->checkMeshset(1, BLOCKSET)) {
438  Range inner_surface;
439  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
440  1, BLOCKSET, 2, inner_surface, true);
441  if (!inner_surface.empty()) {
442  Range inner_surface_verts;
443  CHKERR moab.get_connectivity(inner_surface, inner_surface_verts, false);
444  CHKERR m_field.getInterface<FieldBlas>()->setField(
445  u0, MBVERTEX, inner_surface_verts, "u");
446  }
447  }
448 
449  // Get skin on the body, i.e. body boundary, and apply homogenous Dirichlet
450  // conditions on that boundary.
451  Range surface;
452  CHKERR moab.get_entities_by_dimension(0, 2, surface, false);
453  Skinner skin(&m_field.get_moab());
454  Range edges;
455  CHKERR skin.find_skin(0, surface, false, edges);
456  Range edges_part;
457  ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
458  CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
459  PSTATUS_NOT, -1, &edges_part);
460  Range edges_verts;
461  CHKERR moab.get_connectivity(edges_part, edges_verts, false);
462  // Since Dirichlet b.c. are essential boundary conditions, remove DOFs from
463  // the problem.
464  CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
465  simple_interface->getProblemName(), "u",
466  unite(edges_verts, edges_part));
467 
468  // Create mass matrix, calculate and assemble
469  CHKERR DMCreateMatrix_MoFEM(dm, data->M);
470  CHKERR MatZeroEntries(data->M);
471  boost::shared_ptr<Ele> vol_mass_ele(new Ele(m_field));
472  vol_mass_ele->getOpPtrVector().push_back(new OpAssembleMass(data));
473  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
474  vol_mass_ele);
475  CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
476  CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
477 
478  // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
479  // M^-1G(t,u)
480  data->ksp = createKSP(m_field.get_comm());
481  CHKERR KSPSetOperators(data->ksp, data->M, data->M);
482  CHKERR KSPSetFromOptions(data->ksp);
483  CHKERR KSPSetUp(data->ksp);
484 
485  // Create and setup TS solver
486  auto ts = createTS(m_field.get_comm());
487  // Use IMEX solver, i.e. implicit/explicit solver
488  CHKERR TSSetType(ts, TSARKIMEX);
489  CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
490 
491  // Add element to calculate lhs of stiff part
492  CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
493  vol_ele_stiff_lhs, null, null);
494  // Add element to calculate rhs of stiff part
495  CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
496  vol_ele_stiff_rhs, null, null);
497  // Add element to calculate rhs of slow (nonlinear) part
498  CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
499  vol_ele_slow_rhs, null, null);
500 
501  // Add monitor to time solver
502  boost::shared_ptr<Monitor> monitor_ptr(new Monitor(dm, post_proc));
503  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
504  monitor_ptr, null, null);
505 
506  // Create solution vector
509  CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
510 
511  // Solve problem
512  double ftime = 1;
513  CHKERR TSSetDM(ts, dm);
514  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
515  CHKERR TSSetSolution(ts, X);
516  CHKERR TSSetFromOptions(ts);
517  CHKERR TSSolve(ts, X);
518  }
519  CATCH_ERRORS;
520 
521  // finish work cleaning memory, getting statistics, etc.
523 
524  return 0;
525 }
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:127
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::MatSetValues
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
Definition: EntitiesFieldData.hpp:1631
main
int main(int argc, char *argv[])
Definition: reaction_diffusion.cpp:296
surface
Definition: surface.py:1
ReactionDiffusionEquation::OpAssembleStiffRhs::commonData
boost::shared_ptr< CommonData > commonData
Definition: reaction_diffusion.cpp:191
MoFEM::ProblemsManager
Problem manager is used to build and partition problems.
Definition: ProblemsManager.hpp:21
ReactionDiffusionEquation::OpAssembleSlowRhs::vecF
VectorDouble vecF
Definition: reaction_diffusion.cpp:141
MoFEM::FEMethod
structure for User Loop Methods on finite elements
Definition: LoopMethods.hpp:369
ReactionDiffusionEquation::CommonData::ksp
SmartPetscObj< KSP > ksp
Linear solver.
Definition: reaction_diffusion.cpp:41
MoFEM::createTS
auto createTS(MPI_Comm comm)
Definition: PetscSmartObj.hpp:249
ReactionDiffusionEquation::OpAssembleStiffLhs::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: reaction_diffusion.cpp:208
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::OpSetHOInvJacToScalarBases< 2 >
Definition: HODataOperators.hpp:78
MoFEM::OpCalculateScalarFieldValuesDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
Definition: UserDataOperators.hpp:273
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
ReactionDiffusionEquation
Definition: reaction_diffusion.cpp:12
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:523
BasicFiniteElements.hpp
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::createKSP
auto createKSP(MPI_Comm comm)
Definition: PetscSmartObj.hpp:261
ReactionDiffusionEquation::OpAssembleSlowRhs
Assemble slow part.
Definition: reaction_diffusion.cpp:102
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEM::VecSetValues
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Definition: EntitiesFieldData.hpp:1576
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateHOJac< 2 >
Definition: HODataOperators.hpp:273
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1294
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1975
MoFEM::EntitiesFieldData::EntData::getFTensor0N
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Definition: EntitiesFieldData.hpp:1489
MoFEM::DMCreateGlobalVector_MoFEM
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1167
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
ReactionDiffusionEquation::OpAssembleSlowRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: reaction_diffusion.cpp:105
ReactionDiffusionEquation::OpAssembleStiffLhs::OpAssembleStiffLhs
OpAssembleStiffLhs(boost::shared_ptr< CommonData > &data)
Definition: reaction_diffusion.cpp:204
ReactionDiffusionEquation::EntData
EntitiesFieldData::EntData EntData
Definition: reaction_diffusion.cpp:16
MoFEM::DMCreateMatrix_MoFEM
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1197
ReactionDiffusionEquation::Monitor::preProcess
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Definition: reaction_diffusion.cpp:274
save_every_nth_step
int save_every_nth_step
Definition: photon_diffusion.cpp:67
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: FaceElementForcesAndSourcesCore.hpp:86
Ele
FaceElementForcesAndSourcesCore Ele
Definition: quad_polynomial_approximation.cpp:13
ReactionDiffusionEquation::OpAssembleStiffLhs
Assemble stiff part tangent.
Definition: reaction_diffusion.cpp:203
MoFEM::DMMoFEMTSSetIJacobian
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: DMMoFEM.cpp:853
ReactionDiffusionEquation::OpAssembleMass::commonData
boost::shared_ptr< CommonData > commonData
Definition: reaction_diffusion.cpp:92
ReactionDiffusionEquation::OpAssembleStiffRhs
Assemble stiff part.
Definition: reaction_diffusion.cpp:151
ReactionDiffusionEquation::u0
const double u0
inital vale on blocksets
Definition: reaction_diffusion.cpp:24
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
a
constexpr double a
Definition: approx_sphere.cpp:30
ReactionDiffusionEquation::OpAssembleMass::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Definition: reaction_diffusion.cpp:52
ReactionDiffusionEquation::OpAssembleMass::OpAssembleMass
OpAssembleMass(boost::shared_ptr< CommonData > &data)
Definition: reaction_diffusion.cpp:48
help
static char help[]
Definition: reaction_diffusion.cpp:10
MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator
default operator for TRI element
Definition: FaceElementForcesAndSourcesCore.hpp:94
ReactionDiffusionEquation::OpAssembleStiffLhs::transMat
MatrixDouble transMat
Definition: reaction_diffusion.cpp:260
convert.type
type
Definition: convert.py:64
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
ReactionDiffusionEquation::Monitor::postProc
boost::shared_ptr< PostProcEle > postProc
Definition: reaction_diffusion.cpp:289
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::Simple::addDomainField
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:264
MoFEM::EntitiesFieldData::EntData::getIndices
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Definition: EntitiesFieldData.hpp:1201
surface.surface
def surface(x, y, z, eta)
Definition: surface.py:3
ReactionDiffusionEquation::OpAssembleStiffRhs::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: reaction_diffusion.cpp:154
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
ReactionDiffusionEquation::OpAssembleStiffRhs::vecF
VectorDouble vecF
Definition: reaction_diffusion.cpp:192
ReactionDiffusionEquation::OpAssembleSlowRhs::OpAssembleSlowRhs
OpAssembleSlowRhs(boost::shared_ptr< CommonData > &data)
Definition: reaction_diffusion.cpp:103
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
ReactionDiffusionEquation::OpAssembleStiffLhs::commonData
boost::shared_ptr< CommonData > commonData
Definition: reaction_diffusion.cpp:259
ReactionDiffusionEquation::OpAssembleStiffRhs::OpAssembleStiffRhs
OpAssembleStiffRhs(boost::shared_ptr< CommonData > &data)
Definition: reaction_diffusion.cpp:152
ReactionDiffusionEquation::Monitor
Monitor solution.
Definition: reaction_diffusion.cpp:269
ReactionDiffusionEquation::Monitor::operator()
MoFEMErrorCode operator()()
function is run for every finite element
Definition: reaction_diffusion.cpp:275
FTensor::Index< 'i', DIM >
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
ReactionDiffusionEquation::CommonData
Common data.
Definition: reaction_diffusion.cpp:34
ReactionDiffusionEquation::Monitor::postProcess
MoFEMErrorCode postProcess()
function is run at the end of loop
Definition: reaction_diffusion.cpp:277
Range
MoFEM::DMMoFEMTSSetRHSFunction
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: DMMoFEM.cpp:882
MoFEM::CoreTmp< 0 >::Initialize
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
MoFEM::EntitiesFieldData::EntData::getFTensor1DiffN
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
Definition: EntitiesFieldData.cpp:526
ReactionDiffusionEquation::Monitor::Monitor
Monitor(SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcEle > &post_proc)
Definition: reaction_diffusion.cpp:271
ReactionDiffusionEquation::OpAssembleSlowRhs::commonData
boost::shared_ptr< CommonData > commonData
Definition: reaction_diffusion.cpp:140
ReactionDiffusionEquation::CommonData::dot_val
VectorDouble dot_val
Rate of values of field "u" at integration points.
Definition: reaction_diffusion.cpp:38
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1102
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
BLOCKSET
@ BLOCKSET
Definition: definitions.h:148
ReactionDiffusionEquation::Monitor::dM
SmartPetscObj< DM > dM
Definition: reaction_diffusion.cpp:288
MoFEM::DMMoFEMTSSetMonitor
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: DMMoFEM.cpp:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::DMMoFEMTSSetIFunction
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: DMMoFEM.cpp:800
MoFEM::MeshsetsManager
Interface for managing meshsets containing materials and boundary conditions.
Definition: MeshsetsManager.hpp:104
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3254
ReactionDiffusionEquation::OpAssembleMass
Assemble mass matrix.
Definition: reaction_diffusion.cpp:47
ReactionDiffusionEquation::r
const double r
rate factor
Definition: reaction_diffusion.cpp:21
ReactionDiffusionEquation::CommonData::M
SmartPetscObj< Mat > M
Mass matrix.
Definition: reaction_diffusion.cpp:40
MoFEM::Types::VectorDouble
UBlasVector< double > VectorDouble
Definition: Types.hpp:68
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
ReactionDiffusionEquation::k
const double k
caring capacity
Definition: reaction_diffusion.cpp:22
MoFEM::SmartPetscObj< Mat >
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEM::MeshsetsManager::checkMeshset
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
Definition: MeshsetsManager.cpp:357
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
ReactionDiffusionEquation::CommonData::val
VectorDouble val
Values of field "u" at integration points.
Definition: reaction_diffusion.cpp:37
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
MoFEM::FieldBlas
Basic algebra on fields.
Definition: FieldBlas.hpp:21
ReactionDiffusionEquation::CommonData::grad
MatrixDouble grad
Gradients of field "u" at integration points.
Definition: reaction_diffusion.cpp:36
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
ReactionDiffusionEquation::OpAssembleMass::transMat
MatrixDouble transMat
Definition: reaction_diffusion.cpp:91
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698