v0.9.0
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 OpCalculateScalarValuesDot("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:433
Problem manager is used to build and partition problems.
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:443
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:256
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:434
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:481
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:1042
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:761
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:496
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:1015
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:790
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.
Calculate inverse of jacobian for face element.
MatrixDouble invJac
Inverse of element jacobian.
MoFEMErrorCode loadFile(const std::string options="PARALLEL=READ_PART;" "PARALLEL_RESOLVE_SHARED_ENTS;" "PARTITION=PARALLEL_PARTITION;", const std::string mesh_file_name="")
Load mesh file.
Definition: Simple.cpp:61
Basic algebra on fields.
Definition: FieldBlas.hpp:34
MoFEMErrorCode postProcess()
function is run at the end of loop
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
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:53
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:150
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:105
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:254
#define CHKERR
Inline error check.
Definition: definitions.h:600
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:904
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:233
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Transform local reference derivatives of shape functions to global derivatives.
ForcesAndSourcesCore::UserDataOperator UserDataOperator
boost::shared_ptr< PostProcFaceOnRefinedMesh > postProc
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:289
auto createKSP
Definition: AuxPETSc.hpp:270
Postprocess on face.
static char help[]
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Data on single entity (This is passed as argument to DataOperator::doWork)
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:47
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:708
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:411
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:437
virtual MPI_Comm & get_comm() const =0
continuous field
Definition: definitions.h:175
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:72
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:249
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:61