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