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