v0.13.1
std_reac_diff.cpp
Go to the documentation of this file.
1#include <stdlib.h>
2#include <BasicFiniteElements.hpp>
3#include <rd_stdOperators.hpp>
4
5using namespace MoFEM;
6using namespace StdRDOperators;
7
8static char help[] = "...\n\n";
9
10// struct RhsU {
11// double operator()(const double u, const double v) const {
12// return c * u * (u - alpha) * (1.0 - u) - u * v;
13// }
14// };
15
16// struct RhsV {
17// double operator()(const double u, const double v) const {
18// return (gma + mu1 * v / (mu2 + u)) * (-v - c * u * (u - b - 1.0));
19// }
20// };
21
22struct RDProblem {
23public:
24 RDProblem(moab::Core &mb_instance, MoFEM::Core &core, const int order,
25 const int n_species)
26 : moab(mb_instance)
27 , m_field(core)
28 , order(order)
29 , nb_species(n_species)
30 , cOmm(m_field.get_comm())
31 , rAnk(m_field.get_comm_rank())
32 {
33 vol_ele_slow_rhs = boost::shared_ptr<Ele>(new Ele(m_field));
34 vol_ele_stiff_rhs = boost::shared_ptr<Ele>(new Ele(m_field));
35 vol_ele_stiff_lhs = boost::shared_ptr<Ele>(new Ele(m_field));
36
37 boundary_ele_rhs = boost::shared_ptr<BoundaryEle>(new BoundaryEle(m_field));
38
39 vol_mass_ele = boost::shared_ptr<Ele>(new Ele(m_field));
40
41 post_proc = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
43
44 data.resize(nb_species);
45 values_ptr.resize(nb_species);
46 grads_ptr.resize(nb_species);
47 dots_ptr.resize(nb_species);
49
50 for (int i = 0; i < nb_species; ++i) {
51 data[i] = boost::shared_ptr<PreviousData>(new PreviousData());
52 grads_ptr[i] = boost::shared_ptr<MatrixDouble>(data[i], &data[i]->grads);
53 values_ptr[i] =
54 boost::shared_ptr<VectorDouble>(data[i], &data[i]->values);
55 dots_ptr[i] =
56 boost::shared_ptr<VectorDouble>(data[i], &data[i]->dot_values);
57 }
58 }
59
60 // RDProblem(const int order) : order(order){}
62
63private:
65 MoFEMErrorCode add_fe(std::string field_name);
66 MoFEMErrorCode set_blockData(std::map<int, BlockData> &block_data_map);
67
68 MoFEMErrorCode set_initial_values(std::string field_name, int block_id,
69 Range &surface, double &init_val);
70
71 MoFEMErrorCode update_slow_rhs(std::string mass_fiedl,
72 boost::shared_ptr<VectorDouble> &mass_ptr);
73
74 MoFEMErrorCode push_slow_rhs(std::string mass_name,
75 boost::shared_ptr<PreviousData> &data);
76
78
80
81 MoFEMErrorCode update_vol_fe(boost::shared_ptr<Ele> &vol_ele,
82 boost::shared_ptr<PreviousData> &data);
84 boost::shared_ptr<VectorDouble> &values_ptr,
85 boost::shared_ptr<MatrixDouble> &grads_ptr,
86 boost::shared_ptr<VectorDouble> &dots_ptr);
87
89 boost::shared_ptr<PreviousData> &data,
90 std::map<int, BlockData> &block_map);
91
93 boost::shared_ptr<VectorDouble> &values_ptr,
94 boost::shared_ptr<MatrixDouble> &grads_ptr);
96 boost::shared_ptr<PreviousData> &data,
97 std::map<int, BlockData> &block_map);
98
100
105
108
110
115
117
118 std::vector<Range> inner_surface;
120
121 double global_error;
122
123 MPI_Comm cOmm;
124 const int rAnk;
125
126 int order;
127 int nb_species;
128
129 std::map<int, BlockData> material_blocks;
130
131 boost::shared_ptr<Ele> vol_ele_slow_rhs;
132 boost::shared_ptr<Ele> vol_ele_stiff_rhs;
133 boost::shared_ptr<Ele> vol_ele_stiff_lhs;
134 boost::shared_ptr<BoundaryEle> boundary_ele_rhs;
135
136 boost::shared_ptr<Ele> vol_mass_ele;
137
138 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc;
139 boost::shared_ptr<Monitor> monitor_ptr;
140
141 std::vector<boost::shared_ptr<PreviousData>> data;
142
143 std::vector<boost::shared_ptr<MatrixDouble>> grads_ptr;
144 std::vector<boost::shared_ptr<VectorDouble>> values_ptr;
145 std::vector<boost::shared_ptr<VectorDouble>> dots_ptr;
146
147 boost::shared_ptr<ForcesAndSourcesCore> null;
148};
149
150const double ramp_t = 1.0;
151const double sml = 0.0;
152const double T = 2.0 * M_PI;
153
154struct ExactFunction {
155 double operator()(const double x, const double y, const double t) const {
156 double g = sin(T * x) * sin(T * y);
157 double val = 0;
158 if (x > -sml) {
159 val = 1.0 * g;
160 } else {
161 val = g;
162 }
163 if (t <= ramp_t) {
164 return val * t;
165 } else {
166 return val * ramp_t;
167 }
168 }
169};
170
171struct ExactFunctionGrad {
172 FTensor::Tensor1<double, 3> operator()(const double x, const double y,
173 const double t) const {
175 double mx = -T * cos(T * x) * sin(T * y);
176 double my = -T * sin(T * x) * cos(T * y);
177 double hx, hy;
178 if (x > -sml) {
179 hx = 1.0 * mx;
180 hy = 1.0 * my;
181 } else {
182 hx = mx;
183 hy = my;
184 }
185 if (t <= ramp_t) {
186 grad(0) = hx * t;
187 grad(1) = hy * t;
188 } else {
189 grad(0) = hx * ramp_t;
190 grad(1) = hy * ramp_t;
191 }
192 grad(2) = 0.0;
193 return grad;
194 }
195};
196
197struct ExactFunctionLap {
198 double operator()(const double x, const double y, const double t) const {
199 double glap = -2.0 * pow(T, 2) * sin(T * x) * sin(T * y);
200 double lap;
201 if (x > -sml) {
202 lap = 1.0 * glap;
203 } else {
204 lap = glap;
205 }
206 if (t <= ramp_t) {
207 return lap * t;
208 } else {
209 return lap * ramp_t;
210 }
211 }
212};
213
214struct ExactFunctionDot {
215 double operator()(const double x, const double y, const double t) const {
216 double gdot = sin(T * x) * sin(T * y);
217 double dot;
218 if (x > -sml) {
219 dot = 1.0 * gdot;
220 } else {
221 dot = gdot;
222 }
223 if (t <= ramp_t) {
224 return dot;
225 } else {
226 return 0;
227 }
228 }
229};
230
231// struct ReactionTerms {
232// void operator()(const double u1, const double u2, const double u3,
233// BlockData &coeffs, VectorDouble &vec){
234// vec.resize(3, false);
235// vec.clear();
236
237// vec[0] = coeffs.rate[0] * u1 * (1.0 - coeffs.coef(0, 0) * u1
238// - coeffs.coef(0, 1) * u2
239// - coeffs.coef(0, 2) * u3);
240
241// vec[1] = coeffs.rate[1] * u1 * (1.0 - coeffs.coef(1, 0) * u1
242// - coeffs.coef(1, 1) * u2
243// - coeffs.coef(1, 2) * u3);
244// vec[2] = coeffs.rate[2] * (u1 * u2 + u3) * (1.0 - coeffs.coef(0, 0) * u1
245// - coeffs.coef(0, 1) * u2
246// - coeffs.coef(0, 2) * u3);
247// }
248
249// };
256}
262 1);
263 // CHKERR simple_interface->addBoundaryField(field_name, H1,
264 // AINSWORTH_LEGENDRE_BASE, 1);
267
269}
270
271MoFEMErrorCode RDProblem::set_blockData(std::map<int, BlockData> &block_map) {
274 string name = it->getName();
275 const int id = it->getMeshsetId();
276 if (name.compare(0, 14, "REGION1") == 0) {
277 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
278 id, BLOCKSET, 2, block_map[id].block_ents, true);
279 block_map[id].B0 = 1e-3;
280 block_map[id].block_id = id;
281 } else if (name.compare(0, 14, "REGION2") == 0) {
282 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
283 id, BLOCKSET, 2, block_map[id].block_ents, true);
284 block_map[id].B0 = 1e-2;
285 block_map[id].block_id = id;
286 } else if (name.compare(0, 14, "REGION3") == 0) {
287 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
288 id, BLOCKSET, 2, block_map[id].block_ents, true);
289 block_map[id].B0 = 5e-2;
290 block_map[id].block_id = id;
291 } else if (name.compare(0, 14, "REGION4") == 0) {
292 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
293 id, BLOCKSET, 2, block_map[id].block_ents, true);
294 block_map[id].B0 = 1e-1;
295 block_map[id].block_id = id;
296 } else if (name.compare(0, 14, "REGION5") == 0) {
297 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
298 id, BLOCKSET, 2, block_map[id].block_ents, true);
299 block_map[id].B0 = 5e-1;
300 block_map[id].block_id = id;
301 }
302 }
304}
305
307 int block_id, Range &surface,
308 double &init_val) {
311 BLOCKSET)) {
312 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
313 block_id, BLOCKSET, 2, surface, true);
314 }
315 if (!surface.empty()) {
316 Range surface_verts;
317 CHKERR moab.get_connectivity(surface, surface_verts, false);
319 init_val, MBVERTEX, surface_verts, field_name);
320 }
321
323}
324
326RDProblem::update_slow_rhs(std::string mass_field,
327 boost::shared_ptr<VectorDouble> &values_ptr) {
329 vol_ele_slow_rhs->getOpPtrVector().push_back(
332}
333
335 boost::shared_ptr<PreviousData> &data) {
337
338 vol_ele_slow_rhs->getOpPtrVector().push_back(
342}
343
344MoFEMErrorCode RDProblem::update_vol_fe(boost::shared_ptr<Ele> &vol_ele,
345 boost::shared_ptr<PreviousData> &data) {
347 auto det_ptr = boost::make_shared<VectorDouble>();
348 auto jac_ptr = boost::make_shared<MatrixDouble>();
349 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
350 vol_ele->getOpPtrVector().push_back(new OpCalculateHOJacForFace(jac_ptr));
351 vol_ele->getOpPtrVector().
352 push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
353 vol_ele->getOpPtrVector().push_back(new OpSetInvJacH1ForFace(inv_jac_ptr));
355}
356
359 boost::shared_ptr<VectorDouble> &values_ptr,
360 boost::shared_ptr<MatrixDouble> &grads_ptr,
361 boost::shared_ptr<VectorDouble> &dots_ptr) {
362
364
365 vol_ele_stiff_rhs->getOpPtrVector().push_back(
367 vol_ele_stiff_rhs->getOpPtrVector().push_back(
369 vol_ele_stiff_rhs->getOpPtrVector().push_back(
372}
374 boost::shared_ptr<PreviousData> &data,
375 std::map<int, BlockData> &block_map) {
377
378 vol_ele_stiff_rhs->getOpPtrVector().push_back(
382}
383
386 boost::shared_ptr<VectorDouble> &values_ptr,
387 boost::shared_ptr<MatrixDouble> &grads_ptr) {
389 vol_ele_stiff_lhs->getOpPtrVector().push_back(
391
392 vol_ele_stiff_lhs->getOpPtrVector().push_back(
395}
396
398 boost::shared_ptr<PreviousData> &data,
399 std::map<int, BlockData> &block_map) {
401
402 vol_ele_stiff_lhs->getOpPtrVector().push_back(
404
406}
407
410 auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
411 vol_ele_slow_rhs->getRuleHook = vol_rule;
412
413 vol_ele_stiff_rhs->getRuleHook = vol_rule;
414
415 vol_ele_stiff_lhs->getRuleHook = vol_rule;
416 boundary_ele_rhs->getRuleHook = vol_rule;
418}
419
422 vol_mass_ele->getOpPtrVector().push_back(
425}
426
431 CHKERR MatAssemblyBegin(mass_matrix, MAT_FINAL_ASSEMBLY);
432 CHKERR MatAssemblyEnd(mass_matrix, MAT_FINAL_ASSEMBLY);
433 // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
434 // M^-1G(t,u)
436 CHKERR KSPSetOperators(mass_Ksp, mass_matrix, mass_matrix);
437 CHKERR KSPSetFromOptions(mass_Ksp);
438 CHKERR KSPSetUp(mass_Ksp);
440}
441
444 CHKERR TSSetType(ts, TSARKIMEX);
445 CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
446
449
452
456}
457
460 post_proc->addFieldValuesPostProc(field_name);
462}
463
469}
470
473 // Create solution vector
476 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
477 // Solve problem
478 double ftime = 1;
479 CHKERR TSSetDM(ts, dm);
480 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
481 CHKERR TSSetSolution(ts, X);
482 CHKERR TSSetFromOptions(ts);
483 CHKERR TSSolve(ts, X);
485}
486
489 global_error = 0;
490 std::vector<std::string> mass_names(nb_species);
491
492 for (int i = 0; i < nb_species; ++i) {
493 mass_names[i] = "MASS" + boost::lexical_cast<std::string>(i + 1);
494 }
496
497 for (int i = 0; i < nb_species; ++i) {
498 CHKERR add_fe(mass_names[i]);
499 }
500
502
503 // for (_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(m_field, BLOCKSET, it)) {
504 // string name = it->getName();
505 // if (name.compare(0, 14, "ESSENTIAL") == 0) {
506 // CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
507 // bdry_ents, true);
508 // }
509 // }
510 Range surface;
511 CHKERR moab.get_entities_by_type(0, MBTRI, surface, false);
512 Skinner skin(&m_field.get_moab());
513 Range edges;
514 CHKERR skin.find_skin(0, surface, false, edges);
515 Range edges_part;
516 ParallelComm *pcomm = ParallelComm::get_pcomm(&moab, MYPCOMM_INDEX);
517 CHKERR pcomm->filter_pstatus(edges, PSTATUS_SHARED | PSTATUS_MULTISHARED,
518 PSTATUS_NOT, -1, &edges_part);
519 Range edges_verts;
520 CHKERR moab.get_connectivity(edges_part, edges_verts, false);
521
523
524 VectorDouble initVals;
525 initVals.resize(3, false);
526 initVals.clear();
527
528 initVals[0] = 1.0;
529 initVals[1] = 3;
530 initVals[2] = 0.0;
531
532 for (int i = 0; i < 1; ++i) {
533 CHKERR set_initial_values(mass_names[i], i + 2, inner_surface[i],
534 initVals[i]);
535 CHKERR update_slow_rhs(mass_names[i], values_ptr[i]);
536 }
537
539 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
540 3, BLOCKSET, 2, stimulation_surface, true);
541 }
542
543 if (nb_species == 1) {
544 vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
545 mass_names[0], data[0], data[0], data[0], material_blocks));
546 } else if (nb_species == 2) {
547 vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
548 mass_names[0], data[0], data[1], data[1], material_blocks));
549 } else if (nb_species == 3) {
550 vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
551 mass_names[0], data[0], data[1], data[2], material_blocks));
552 }
553
554 for (int i = 0; i < nb_species; ++i) {
555 CHKERR push_slow_rhs(mass_names[i], data[i]);
556 // boundary_ele_rhs->getOpPtrVector().push_back(
557 // new OpAssembleNaturalBCRhs(mass_names[i], bdry_ents));
558 }
559
561
562 // cout << "essential : " << essential_bdry_ents.empty() << endl;
563
564 auto solve_for_g = [&]() {
566 if (vol_ele_slow_rhs->vecAssembleSwitch) {
567 CHKERR VecGhostUpdateBegin(vol_ele_slow_rhs->ts_F, ADD_VALUES,
568 SCATTER_REVERSE);
569 CHKERR VecGhostUpdateEnd(vol_ele_slow_rhs->ts_F, ADD_VALUES,
570 SCATTER_REVERSE);
571 CHKERR VecAssemblyBegin(vol_ele_slow_rhs->ts_F);
572 CHKERR VecAssemblyEnd(vol_ele_slow_rhs->ts_F);
573 *vol_ele_slow_rhs->vecAssembleSwitch = false;
574 }
575 CHKERR KSPSolve(mass_Ksp, vol_ele_slow_rhs->ts_F, vol_ele_slow_rhs->ts_F);
577 };
578 // Add hook to the element to calculate g.
579 bdry_ents = unite(edges_verts, edges_part);
580
581 // CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
582 // simple_interface->getProblemName(), mass_names[0], bdry_ents);
583
585 CHKERR MatZeroEntries(mass_matrix);
586
587 for (int i = 0; i < nb_species; ++i) {
588 CHKERR push_mass_ele(mass_names[i]);
589 }
590
592
594
595 for (int i = 0; i < nb_species; ++i) {
597 dots_ptr[i]);
599 }
600
602
603 for (int i = 0; i < nb_species; ++i) {
605 CHKERR push_stiff_lhs(mass_names[i], data[i],
606 material_blocks); // nb_species times
607 }
608 // vol_ele_stiff_lhs->getOpPtrVector().push_back(
609 // new OpError(ExactFunction(), ExactFunctionLap(), ExactFunctionGrad(),
610 // data[0], material_blocks, global_error));
611
613
615
616 vol_ele_slow_rhs->postProcessHook = solve_for_g;
617
619
620 post_proc->generateReferenceElementMesh(); // only once
621
622 for (int i = 0; i < nb_species; ++i) {
623 CHKERR post_proc_fields(mass_names[i]);
624 }
625 post_proc->addFieldValuesPostProc("ERROR");
626
627 monitor_ptr = boost::shared_ptr<Monitor>(
628 new Monitor(cOmm, rAnk, dm, post_proc, global_error)); // nb_species times
629 CHKERR output_result(); // only once
630
631 CHKERR solve(); // only once
632
634}
635
636int main(int argc, char *argv[]) {
637 const char param_file[] = "param_file.petsc";
639 try {
640 moab::Core mb_instance;
641 moab::Interface &moab = mb_instance;
642 MoFEM::Core core(moab);
643 DMType dm_name = "DMMOFEM";
644 CHKERR DMRegister_MoFEM(dm_name);
645
646 int order = 1;
647 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
648 int nb_species = 1;
649 RDProblem reac_diff_problem(mb_instance, core, order, nb_species);
650 CHKERR reac_diff_problem.run_analysis();
651 }
654 return 0;
655}
static Index< 'p', 3 > p
std::string param_file
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ 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
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
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)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpCalculateScalarFieldValuesDot OpCalculateScalarValuesDot
auto createTS(MPI_Comm comm)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
constexpr double t
plate stiffness
Definition: plate.cpp:76
constexpr auto field_name
FaceElementForcesAndSourcesCore Ele
constexpr double g
int main(int argc, char *argv[])
static char help[]
const double sml
const double ramp_t
const double T
double operator()(const double x, const double y, const double t) const
Exact gradient.
FTensor::Tensor1< double, 3 > operator()(const double x, const double y, const double t) const
double operator()(const double x, const double y, const double t) const
double operator()(const double x, const double y, const double t) const
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
Deprecated interface functions.
Basic algebra on fields.
Definition: FieldBlas.hpp:31
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.
Simple interface for fast problem set-up.
Definition: Simple.hpp:35
MoFEMErrorCode addDataField(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:434
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 getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:321
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Postprocess on face.
std::vector< Range > inner_surface
MoFEMErrorCode update_slow_rhs(std::string mass_fiedl, boost::shared_ptr< VectorDouble > &mass_ptr)
moab::Interface & moab
boost::shared_ptr< FaceEle > vol_ele_stiff_lhs
std::map< int, BlockData > material_blocks
SmartPetscObj< TS > ts
RDProblem(moab::Core &mb_instance, MoFEM::Core &core, const int order, const int n_species)
MoFEMErrorCode update_stiff_rhs()
SmartPetscObj< DM > dm
boost::shared_ptr< Ele > vol_ele_stiff_lhs
boost::shared_ptr< Ele > vol_mass_ele
Range stimulation_surface
MoFEMErrorCode set_fe_in_loop()
boost::shared_ptr< Monitor > monitor_ptr
MoFEMErrorCode add_fe()
SmartPetscObj< Mat > mass_matrix
MoFEMErrorCode push_slow_rhs()
MoFEMErrorCode update_vol_fe(boost::shared_ptr< FaceEle > &vol_ele, boost::shared_ptr< PreviousData > &data)
boost::shared_ptr< Ele > vol_ele_slow_rhs
MoFEMErrorCode post_proc_fields()
boost::shared_ptr< ForcesAndSourcesCore > null
boost::shared_ptr< BoundaryEle > boundary_ele_rhs
MoFEM::Interface & m_field
Simple * simple_interface
MoFEMErrorCode run_analysis()
MPI_Comm cOmm
boost::shared_ptr< Ele > vol_ele_stiff_rhs
double global_error
MoFEMErrorCode set_blockData(std::map< int, BlockData > &block_data_map)
SmartPetscObj< KSP > mass_Ksp
std::vector< boost::shared_ptr< VectorDouble > > values_ptr
const int rAnk
MoFEMErrorCode resolve_slow_rhs()
MoFEMErrorCode push_stiff_lhs()
MoFEMErrorCode push_stiff_rhs()
MoFEMErrorCode solve()
std::vector< boost::shared_ptr< VectorDouble > > dots_ptr
MoFEMErrorCode set_initial_values(std::string field_name, int block_id, Range &surface, double &init_val)
MoFEMErrorCode push_mass_ele(std::string field_name)
MoFEMErrorCode setup_system()
MoFEMErrorCode output_result()
std::vector< boost::shared_ptr< PreviousData > > data
boost::shared_ptr< FaceEle > vol_ele_slow_rhs
MoFEMErrorCode update_stiff_lhs()
std::vector< boost::shared_ptr< MatrixDouble > > grads_ptr
MoFEMErrorCode set_integration_rule()
boost::shared_ptr< PostProcFaceOnRefinedMesh > post_proc
boost::shared_ptr< FaceEle > vol_ele_stiff_rhs