v0.14.0
Loading...
Searching...
No Matches
reaction_diffusion.cpp
Go to the documentation of this file.
1/**
2 * \file reaction_diffusion.cpp
3 * \example reaction_diffusion.cpp
4 *
5 **/
6
8using namespace MoFEM;
9
10static char help[] = "...\n\n";
11
13
17
19
20const double D = 2e-3; ///< diffusivity
21const double r = 1; ///< rate factor
22const double k = 1; ///< caring capacity
23
24const double u0 = 0.1; ///< inital vale on blocksets
25
26const int save_every_nth_step = 1;
27
28/**
29 * @brief Common data
30 *
31 * Common data are used to keep and pass data between elements
32 *
33 */
34struct CommonData {
35
36 MatrixDouble grad; ///< Gradients of field "u" at integration points
37 VectorDouble val; ///< Values of field "u" at integration points
38 VectorDouble dot_val; ///< Rate of values of field "u" at integration points
39
40 SmartPetscObj<Mat> M; ///< Mass matrix
41 SmartPetscObj<KSP> ksp; ///< Linear solver
42};
43
44/**
45 * @brief Assemble mass matrix
46 */
48 OpAssembleMass(boost::shared_ptr<CommonData> &data)
49 : OpEle("u", "u", OpEle::OPROWCOL), commonData(data) {
50 sYmm = true;
51 }
52 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
53 EntityType col_type, EntData &row_data,
54 EntData &col_data) {
56 const int nb_row_dofs = row_data.getIndices().size();
57 const int nb_col_dofs = col_data.getIndices().size();
58 if (nb_row_dofs && nb_col_dofs) {
59 const int nb_integration_pts = getGaussPts().size2();
60 mat.resize(nb_row_dofs, nb_col_dofs, false);
61 mat.clear();
62 auto t_row_base = row_data.getFTensor0N();
64 const double vol = getMeasure();
65 for (int gg = 0; gg != nb_integration_pts; ++gg) {
66 const double a = t_w * vol;
67 for (int rr = 0; rr != nb_row_dofs; ++rr) {
68 auto t_col_base = col_data.getFTensor0N(gg, 0);
69 for (int cc = 0; cc != nb_col_dofs; ++cc) {
70 mat(rr, cc) += a * t_row_base * t_col_base;
71 ++t_col_base;
72 }
73 ++t_row_base;
74 }
75 ++t_w;
76 }
77
78 CHKERR MatSetValues(commonData->M, row_data, col_data, &mat(0, 0),
79 ADD_VALUES);
80 if (row_side != col_side || row_type != col_type) {
81 transMat.resize(nb_col_dofs, nb_row_dofs, false);
82 noalias(transMat) = trans(mat);
83 CHKERR MatSetValues(commonData->M, col_data, row_data, &transMat(0, 0),
84 ADD_VALUES);
85 }
86 }
88 }
89
90private:
92 boost::shared_ptr<CommonData> commonData;
93};
94
95/**
96 * @brief Assemble slow part
97 *
98 * Solve problem \f$ F(t,u,\dot{u}) = G(t,u) \f$ where here the right hand side
99 * \f$ G(t,u) \f$ is implemented.
100 *
101 */
103 OpAssembleSlowRhs(boost::shared_ptr<CommonData> &data)
104 : OpEle("u", OpEle::OPROW), commonData(data) {}
105 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
107 const int nb_dofs = data.getIndices().size();
108 if (nb_dofs) {
109 vecF.resize(nb_dofs, false);
110 vecF.clear();
111
112 const int nb_integration_pts = getGaussPts().size2();
113 auto t_val = getFTensor0FromVec(commonData->val);
114 auto t_base = data.getFTensor0N();
115 auto t_w = getFTensor0IntegrationWeight();
116
117 const double vol = getMeasure();
118 for (int gg = 0; gg != nb_integration_pts; ++gg) {
119 const double a = vol * t_w;
120 const double f = a * r * t_val * (1 - t_val / k);
121 for (int rr = 0; rr != nb_dofs; ++rr) {
122 const double b = f * t_base;
123 vecF[rr] += b;
124 ++t_base;
125 }
126
127 ++t_val;
128 ++t_w;
129 }
130
131 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
132 PETSC_TRUE);
133 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
134 ADD_VALUES);
135 }
137 }
138
139private:
140 boost::shared_ptr<CommonData> commonData;
142};
143
144/**
145 * @brief Assemble stiff part
146 *
147 * Solve problem \f$ F(t,u,\dot{u}) = G(t,u) \f$ where here the right hand side
148 * \f$ F(t,u,\dot{u}) \f$ is implemented.
149 *
150 */
151template <int DIM> struct OpAssembleStiffRhs : OpEle {
152 OpAssembleStiffRhs(boost::shared_ptr<CommonData> &data)
153 : OpEle("u", OpEle::OPROW), commonData(data) {}
154 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
156 const int nb_dofs = data.getIndices().size();
157 if (nb_dofs) {
158 vecF.resize(nb_dofs, false);
159 vecF.clear();
160
161 const int nb_integration_pts = getGaussPts().size2();
162 auto t_dot_val = getFTensor0FromVec(commonData->dot_val);
163 auto t_grad = getFTensor1FromMat<DIM>(commonData->grad);
164 auto t_base = data.getFTensor0N();
165 auto t_diff_base = data.getFTensor1DiffN<DIM>();
166 auto t_w = getFTensor0IntegrationWeight();
168
169 const double vol = getMeasure();
170 for (int gg = 0; gg != nb_integration_pts; ++gg) {
171 const double a = vol * t_w;
172 for (int rr = 0; rr != nb_dofs; ++rr) {
173 vecF[rr] += a * (t_base * t_dot_val + D * t_diff_base(i) * t_grad(i));
174 ++t_diff_base;
175 ++t_base;
176 }
177 ++t_dot_val;
178 ++t_grad;
179 ++t_w;
180 }
181
182 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
183 PETSC_TRUE);
184 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
185 ADD_VALUES);
186 }
188 }
189
190private:
191 boost::shared_ptr<CommonData> commonData;
193};
194
195/**
196 * @brief Assemble stiff part tangent
197 *
198 * Solve problem \f$ F(t,u,\dot{u}) = G(t,u) \f$ where here the right hand side
199 * \f$ \frac{\textrm{d} F}{\textrm{d} u^n} = a F_{\dot{u}}(t,u,\textrm{u}) +
200 * F_{u}(t,u,\textrm{u}) \f$ is implemented.
201 *
202 */
203template <int DIM> struct OpAssembleStiffLhs : OpEle {
204 OpAssembleStiffLhs(boost::shared_ptr<CommonData> &data)
205 : OpEle("u", "u", OpEle::OPROWCOL), commonData(data) {
206 sYmm = true;
207 }
208 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
209 EntityType col_type, EntData &row_data,
210 EntData &col_data) {
212 const int nb_row_dofs = row_data.getIndices().size();
213 const int nb_col_dofs = col_data.getIndices().size();
214 if (nb_row_dofs && nb_col_dofs) {
215 mat.resize(nb_row_dofs, nb_col_dofs, false);
216 mat.clear();
217
218 const int nb_integration_pts = getGaussPts().size2();
219 auto t_row_base = row_data.getFTensor0N();
220 auto t_row_diff_base = row_data.getFTensor1DiffN<DIM>();
221 auto t_w = getFTensor0IntegrationWeight();
223
224 const double ts_a = getFEMethod()->ts_a;
225 const double vol = getMeasure();
226 for (int gg = 0; gg != nb_integration_pts; ++gg) {
227 const double a = vol * t_w;
228 for (int rr = 0; rr != nb_row_dofs; ++rr) {
229
230 auto t_col_base = col_data.getFTensor0N(gg, 0);
231 auto t_col_diff_base = col_data.getFTensor1DiffN<DIM>(gg, 0);
232
233 for (int cc = 0; cc != nb_col_dofs; ++cc) {
234 mat(rr, cc) += a * (t_row_base * t_col_base * ts_a +
235 D * t_row_diff_base(i) * t_col_diff_base(i));
236 ++t_col_base;
237 ++t_col_diff_base;
238 }
239
240 ++t_row_base;
241 ++t_row_diff_base;
242 }
243 ++t_w;
244 }
245
246 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
247 ADD_VALUES);
248 if (row_side != col_side || row_type != col_type) {
249 transMat.resize(nb_col_dofs, nb_row_dofs, false);
250 noalias(transMat) = trans(mat);
251 CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
252 &transMat(0, 0), ADD_VALUES);
253 }
254 }
256 }
257
258private:
259 boost::shared_ptr<CommonData> commonData;
261};
262
263/**
264 * @brief Monitor solution
265 *
266 * This functions is called by TS solver at the end of each step. It is used
267 * to output results to the hard drive.
268 */
269struct Monitor : public FEMethod {
270
271 Monitor(SmartPetscObj<DM> &dm, boost::shared_ptr<PostProcEle> &post_proc)
272 : dM(dm), postProc(post_proc){};
273
276
279 if (ts_step % save_every_nth_step == 0) {
281 CHKERR postProc->writeFile(
282 "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
283 }
285 }
286
287private:
289 boost::shared_ptr<PostProcEle> postProc;
290};
291
292}; // namespace ReactionDiffusionEquation
293
294using namespace ReactionDiffusionEquation;
295
296int main(int argc, char *argv[]) {
297
298 // initialize petsc
299 const char param_file[] = "param_file.petsc";
301
302 try {
303
304 // Create moab and mofem instances
305 moab::Core mb_instance;
306 moab::Interface &moab = mb_instance;
307 MoFEM::Core core(moab);
308 MoFEM::Interface &m_field = core;
309
310 // Register DM Manager
311 DMType dm_name = "DMMOFEM";
312 CHKERR DMRegister_MoFEM(dm_name);
313
314 // Simple interface
315 Simple *simple_interface;
316 CHKERR m_field.getInterface(simple_interface);
317 CHKERR simple_interface->getOptions();
318 CHKERR simple_interface->loadFile();
319
320 int order = 4; ///< approximation order
321 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
322
323 // add fields
324 CHKERR simple_interface->addDomainField("u", H1, AINSWORTH_LEGENDRE_BASE,
325 1);
326 // set fields order
327 CHKERR simple_interface->setFieldOrder("u", order);
328 // setup problem
329 CHKERR simple_interface->setUp();
330
331 // Create common data structure
332 boost::shared_ptr<CommonData> data(new CommonData());
333 /// Alias pointers to data in common data structure
334 auto val_ptr = boost::shared_ptr<VectorDouble>(data, &data->val);
335 auto dot_val_ptr = boost::shared_ptr<VectorDouble>(data, &data->dot_val);
336 auto grad_ptr = boost::shared_ptr<MatrixDouble>(data, &data->grad);
337
338 // Create finite element instances to integrate the right-hand side of slow
339 // and stiff vector, and the tangent left-hand side for stiff part.
340 boost::shared_ptr<Ele> vol_ele_slow_rhs(new Ele(m_field));
341 boost::shared_ptr<Ele> vol_ele_stiff_rhs(new Ele(m_field));
342 boost::shared_ptr<Ele> vol_ele_stiff_lhs(new Ele(m_field));
343
344 // Push operators to integrate the slow right-hand side vector
345 vol_ele_slow_rhs->getOpPtrVector().push_back(
346 new OpCalculateScalarFieldValues("u", val_ptr));
347 vol_ele_slow_rhs->getOpPtrVector().push_back(new OpAssembleSlowRhs(data));
348
349 // PETSc IMAX and Explicit solver demans that g = M^-1 G is provided. So
350 // when the slow right-hand side vector (G) is assembled is solved for g
351 // vector.
352 auto solve_for_g = [&]() {
354 if (*(vol_ele_slow_rhs->vecAssembleSwitch)) {
355 CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
356 SCATTER_REVERSE);
357 CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
358 SCATTER_REVERSE);
359 CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
360 CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
361 *vol_ele_slow_rhs->vecAssembleSwitch = false;
362 }
363 CHKERR KSPSolve(data->ksp, vol_ele_slow_rhs->ts_F,
364 vol_ele_slow_rhs->ts_F);
366 };
367 // Add hook to the element to calculate g.
368 vol_ele_slow_rhs->postProcessHook = solve_for_g;
369
370 auto det_ptr = boost::make_shared<VectorDouble>();
371 auto jac_ptr = boost::make_shared<MatrixDouble>();
372 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
373 // Add operators to calculate the stiff right-hand side
374 vol_ele_stiff_rhs->getOpPtrVector().push_back(
375 new OpCalculateHOJac<2>(jac_ptr));
376 vol_ele_stiff_rhs->getOpPtrVector().push_back(
377 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
378 vol_ele_stiff_rhs->getOpPtrVector().push_back(
379 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
380 vol_ele_stiff_rhs->getOpPtrVector().push_back(
381 new OpCalculateScalarFieldValuesDot("u", dot_val_ptr));
382 vol_ele_stiff_rhs->getOpPtrVector().push_back(
383 new OpCalculateScalarFieldGradient<2>("u", grad_ptr));
384 vol_ele_stiff_rhs->getOpPtrVector().push_back(
385 new OpAssembleStiffRhs<2>(data));
386
387 // Add operators to calculate the stiff left-hand side
388 vol_ele_stiff_lhs->getOpPtrVector().push_back(
389 new OpCalculateHOJac<2>(jac_ptr));
390 vol_ele_stiff_lhs->getOpPtrVector().push_back(
391 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
392 vol_ele_stiff_lhs->getOpPtrVector().push_back(
393 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
394 vol_ele_stiff_lhs->getOpPtrVector().push_back(
395 new OpAssembleStiffLhs<2>(data));
396
397 // Set integration rules
398 auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
399 vol_ele_slow_rhs->getRuleHook = vol_rule;
400 vol_ele_stiff_rhs->getRuleHook = vol_rule;
401 vol_ele_stiff_lhs->getRuleHook = vol_rule;
402
403 // Crate element for post-processing
404
405 auto post_proc = boost::make_shared<PostProcEle>(m_field);
406 boost::shared_ptr<ForcesAndSourcesCore> null;
407
409
410 auto u_ptr = boost::make_shared<VectorDouble>();
411 post_proc->getOpPtrVector().push_back(
412 new OpCalculateScalarFieldValues("u", u_ptr));
413 post_proc->getOpPtrVector().push_back(
414
415 new OpPPMap(
416
417 post_proc->getPostProcMesh(), post_proc->getMapGaussPts(),
418
419 {{"u", u_ptr}},
420
421 {},
422
423 {},
424
425 {}
426
427 )
428
429 );
430
431 // Get PETSc discrete manager
432 auto dm = simple_interface->getDM();
433
434 // Get surface entities form blockset, set initial values in those
435 // blocksets. To keep it simple is assumed that inital values are on
436 // blockset 1
437 if (m_field.getInterface<MeshsetsManager>()->checkMeshset(1, BLOCKSET)) {
438 Range inner_surface;
439 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
440 1, BLOCKSET, 2, inner_surface, true);
441 if (!inner_surface.empty()) {
442 Range inner_surface_verts;
443 CHKERR moab.get_connectivity(inner_surface, inner_surface_verts, false);
444 CHKERR m_field.getInterface<FieldBlas>()->setField(
445 u0, MBVERTEX, inner_surface_verts, "u");
446 }
447 }
448
449 // Get skin on the body, i.e. body boundary, and apply homogenous Dirichlet
450 // conditions on that boundary.
452 CHKERR moab.get_entities_by_dimension(0, 2, surface, false);
453 Skinner skin(&m_field.get_moab());
454 Range edges;
455 CHKERR skin.find_skin(0, surface, false, edges);
456 Range edges_part;
457 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
458 CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
459 PSTATUS_NOT, -1, &edges_part);
460 Range edges_verts;
461 CHKERR moab.get_connectivity(edges_part, edges_verts, false);
462 // Since Dirichlet b.c. are essential boundary conditions, remove DOFs from
463 // the problem.
464 CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
465 simple_interface->getProblemName(), "u",
466 unite(edges_verts, edges_part));
467
468 // Create mass matrix, calculate and assemble
469 CHKERR DMCreateMatrix_MoFEM(dm, data->M);
470 CHKERR MatZeroEntries(data->M);
471 boost::shared_ptr<Ele> vol_mass_ele(new Ele(m_field));
472 vol_mass_ele->getOpPtrVector().push_back(new OpAssembleMass(data));
473 CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
474 vol_mass_ele);
475 CHKERR MatAssemblyBegin(data->M, MAT_FINAL_ASSEMBLY);
476 CHKERR MatAssemblyEnd(data->M, MAT_FINAL_ASSEMBLY);
477
478 // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
479 // M^-1G(t,u)
480 data->ksp = createKSP(m_field.get_comm());
481 CHKERR KSPSetOperators(data->ksp, data->M, data->M);
482 CHKERR KSPSetFromOptions(data->ksp);
483 CHKERR KSPSetUp(data->ksp);
484
485 // Create and setup TS solver
486 auto ts = createTS(m_field.get_comm());
487 // Use IMEX solver, i.e. implicit/explicit solver
488 CHKERR TSSetType(ts, TSARKIMEX);
489 CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
490
491 // Add element to calculate lhs of stiff part
492 CHKERR DMMoFEMTSSetIJacobian(dm, simple_interface->getDomainFEName(),
493 vol_ele_stiff_lhs, null, null);
494 // Add element to calculate rhs of stiff part
495 CHKERR DMMoFEMTSSetIFunction(dm, simple_interface->getDomainFEName(),
496 vol_ele_stiff_rhs, null, null);
497 // Add element to calculate rhs of slow (nonlinear) part
498 CHKERR DMMoFEMTSSetRHSFunction(dm, simple_interface->getDomainFEName(),
499 vol_ele_slow_rhs, null, null);
500
501 // Add monitor to time solver
502 boost::shared_ptr<Monitor> monitor_ptr(new Monitor(dm, post_proc));
503 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple_interface->getDomainFEName(),
504 monitor_ptr, null, null);
505
506 // Create solution vector
509 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
510
511 // Solve problem
512 double ftime = 1;
513 CHKERR TSSetDM(ts, dm);
514 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
515 CHKERR TSSetSolution(ts, X);
516 CHKERR TSSetFromOptions(ts);
517 CHKERR TSSolve(ts, X);
518 }
520
521 // finish work cleaning memory, getting statistics, etc.
523
524 return 0;
525}
static Index< 'p', 3 > p
std::string param_file
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ H1
continuous field
Definition: definitions.h:85
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr int order
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: DMMoFEM.cpp:786
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1183
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:572
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: DMMoFEM.cpp:839
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMoFEM.cpp:1153
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:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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: DMMoFEM.cpp:1042
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: DMMoFEM.cpp:868
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createTS(MPI_Comm comm)
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
const double D
diffusivity
const double r
rate factor
const double u0
inital vale on blocksets
const double k
caring capacity
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FaceElementForcesAndSourcesCore Ele
static char help[]
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
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:21
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ 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.
Post post-proc data at points from hash maps.
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:27
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:264
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
intrusive_ptr for managing petsc objects
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step number
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
[Push operators to pipeline]
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< PostProcEle > 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< PostProcEle > &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.