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