v0.13.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>
22using namespace MoFEM;
23
24static char help[] = "...\n\n";
25
27
31
32const double D = 2e-3; ///< diffusivity
33const double r = 1; ///< rate factor
34const double k = 1; ///< caring capacity
35
36const double u0 = 0.1; ///< inital vale on blocksets
37
38const 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 */
46struct 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();
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
102private:
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) {}
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
151private:
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 */
163template <int DIM> struct OpAssembleStiffRhs : OpEle {
164 OpAssembleStiffRhs(boost::shared_ptr<CommonData> &data)
165 : OpEle("u", OpEle::OPROW), commonData(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
202private:
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 */
215template <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
270private:
271 boost::shared_ptr<CommonData> commonData;
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 */
281struct Monitor : public FEMethod {
282
284 boost::shared_ptr<PostProcFaceOnRefinedMesh> &post_proc)
285 : dM(dm), postProc(post_proc){};
286
289
292 if (ts_step % save_every_nth_step == 0) {
294 CHKERR postProc->writeFile(
295 "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
296 }
298 }
299
300private:
302 boost::shared_ptr<PostProcFaceOnRefinedMesh> postProc;
303};
304
305}; // namespace ReactionDiffusionEquation
306
307using namespace ReactionDiffusionEquation;
308
309int main(int argc, char *argv[]) {
310
311 // initialize petsc
312 const char param_file[] = "param_file.petsc";
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 OpCalculateHOJac<2>(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 OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
393 vol_ele_stiff_rhs->getOpPtrVector().push_back(new OpSetHOWeights(det_ptr));
394 vol_ele_stiff_rhs->getOpPtrVector().push_back(
395 new OpCalculateScalarFieldValuesDot("u", dot_val_ptr));
396 vol_ele_stiff_rhs->getOpPtrVector().push_back(
397 new OpCalculateScalarFieldGradient<2>("u", grad_ptr));
398 vol_ele_stiff_rhs->getOpPtrVector().push_back(
399 new OpAssembleStiffRhs<2>(data));
400
401 // Add operators to calculate the stiff left-hand side
402 vol_ele_stiff_lhs->getOpPtrVector().push_back(
403 new OpCalculateHOJac<2>(jac_ptr));
404 vol_ele_stiff_lhs->getOpPtrVector().push_back(
405 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
406 vol_ele_stiff_lhs->getOpPtrVector().push_back(
407 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
408 vol_ele_stiff_lhs->getOpPtrVector().push_back(new OpSetHOWeights(det_ptr));
409 vol_ele_stiff_lhs->getOpPtrVector().push_back(
410 new OpAssembleStiffLhs<2>(data));
411
412 // Set integration rules
413 auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
414 vol_ele_slow_rhs->getRuleHook = vol_rule;
415 vol_ele_stiff_rhs->getRuleHook = vol_rule;
416 vol_ele_stiff_lhs->getRuleHook = vol_rule;
417
418 // Crate element for post-processing
419 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc =
420 boost::shared_ptr<PostProcFaceOnRefinedMesh>(
421 new PostProcFaceOnRefinedMesh(m_field));
422 boost::shared_ptr<ForcesAndSourcesCore> null;
423 // Genarte post-processing mesh
424 post_proc->generateReferenceElementMesh();
425 // Postprocess only field values
426 post_proc->addFieldValuesPostProc("u");
427
428 // Get PETSc discrete manager
429 auto dm = simple_interface->getDM();
430
431 // Get surface entities form blockset, set initial values in those
432 // blocksets. To keep it simple is assumed that inital values are on
433 // blockset 1
435 Range inner_surface;
436 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
437 1, BLOCKSET, 2, inner_surface, true);
438 if (!inner_surface.empty()) {
439 Range inner_surface_verts;
440 CHKERR moab.get_connectivity(inner_surface, inner_surface_verts, false);
441 CHKERR m_field.getInterface<FieldBlas>()->setField(
442 u0, MBVERTEX, inner_surface_verts, "u");
443 }
444 }
445
446 // Get skin on the body, i.e. body boundary, and apply homogenous Dirichlet
447 // conditions on that boundary.
448 Range surface;
449 CHKERR moab.get_entities_by_dimension(0, 2, surface, false);
450 Skinner skin(&m_field.get_moab());
451 Range edges;
452 CHKERR skin.find_skin(0, surface, false, edges);
453 Range edges_part;
454 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
455 CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
456 PSTATUS_NOT, -1, &edges_part);
457 Range edges_verts;
458 CHKERR moab.get_connectivity(edges_part, edges_verts, false);
459 // Since Dirichlet b.c. are essential boundary conditions, remove DOFs from
460 // the problem.
461 CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
462 simple_interface->getProblemName(), "u",
463 unite(edges_verts, edges_part));
464
465 // Create mass matrix, calculate and assemble
466 CHKERR DMCreateMatrix_MoFEM(dm, data->M);
467 CHKERR MatZeroEntries(data->M);
468 boost::shared_ptr<Ele> vol_mass_ele(new Ele(m_field));
469 vol_mass_ele->getOpPtrVector().push_back(new OpAssembleMass(data));
470 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
471 vol_mass_ele);
472 CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
473 CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
474
475 // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
476 // M^-1G(t,u)
477 data->ksp = createKSP(m_field.get_comm());
478 CHKERR KSPSetOperators(data->ksp, data->M, data->M);
479 CHKERR KSPSetFromOptions(data->ksp);
480 CHKERR KSPSetUp(data->ksp);
481
482 // Create and setup TS solver
483 auto ts = createTS(m_field.get_comm());
484 // Use IMEX solver, i.e. implicit/explicit solver
485 CHKERR TSSetType(ts, TSARKIMEX);
486 CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
487
488 // Add element to calculate lhs of stiff part
489 CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
490 vol_ele_stiff_lhs, null, null);
491 // Add element to calculate rhs of stiff part
492 CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
493 vol_ele_stiff_rhs, null, null);
494 // Add element to calculate rhs of slow (nonlinear) part
495 CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
496 vol_ele_slow_rhs, null, null);
497
498 // Add monitor to time solver
499 boost::shared_ptr<Monitor> monitor_ptr(new Monitor(dm, post_proc));
500 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
501 monitor_ptr, null, null);
502
503 // Create solution vector
506 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
507
508 // Solve problem
509 double ftime = 1;
510 CHKERR TSSetDM(ts, dm);
511 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
512 CHKERR TSSetSolution(ts, X);
513 CHKERR TSSetFromOptions(ts);
514 CHKERR TSSolve(ts, X);
515 }
517
518 // finish work cleaning memory, getting statistics, etc.
520
521 return 0;
522}
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:759
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:482
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1156
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:545
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:812
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1126
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
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
auto createKSP(MPI_Comm comm)
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:1015
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:841
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createTS(MPI_Comm comm)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
OpPlasticTools::CommonData CommonData
EntitiesFieldData::EntData EntData
const double D
diffusivity
const double r
rate factor
const double u0
inital vale on blocksets
const double k
caring capacity
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
bool sYmm
If true assume that matrix is symmetric structure.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
structure for User Loop Methods on finite elements
Basic algebra on fields.
Definition: FieldBlas.hpp:31
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
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.
Set inverse jacobian to base functions.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:35
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:378
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:294
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:810
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:308
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:593
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:753
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:342
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:321
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step
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.