v0.13.1
mixed_reac_diff.cpp
Go to the documentation of this file.
1#include <stdlib.h>
2#include <BasicFiniteElements.hpp>
3#include <RDOperators.hpp>
4
5using namespace MoFEM;
6using namespace ReactionDiffusion;
7
8static char help[] = "...\n\n";
9
10// #define M_PI 3.14159265358979323846 /* pi */
11
12struct RDProblem {
13public:
14 RDProblem(MoFEM::Core &core, const int order)
15 : m_field(core)
16 , order(order)
17 , cOmm(m_field.get_comm())
18 , rAnk(m_field.get_comm_rank()) {
19 vol_ele_slow_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
21 boost::shared_ptr<BoundaryEle>(new BoundaryEle(m_field));
22 vol_ele_stiff_rhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
23 vol_ele_stiff_lhs = boost::shared_ptr<FaceEle>(new FaceEle(m_field));
24 post_proc = boost::shared_ptr<PostProcFaceOnRefinedMesh>(
26
27 data1 = boost::shared_ptr<PreviousData>(new PreviousData());
28 data2 = boost::shared_ptr<PreviousData>(new PreviousData());
29 data3 = boost::shared_ptr<PreviousData>(new PreviousData());
30
32 boost::shared_ptr<MatrixDouble>(data1, &data1->flux_values);
33
34 flux_divs_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->flux_divs);
35
37 boost::shared_ptr<VectorDouble>(data1, &data1->mass_values);
38
39 mass_dots_ptr1 = boost::shared_ptr<VectorDouble>(data1, &data1->mass_dots);
40
42 boost::shared_ptr<MatrixDouble>(data2, &data2->flux_values);
43 flux_divs_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->flux_divs);
44
46 boost::shared_ptr<VectorDouble>(data2, &data2->mass_values);
47 mass_dots_ptr2 = boost::shared_ptr<VectorDouble>(data2, &data2->mass_dots);
48
50 boost::shared_ptr<MatrixDouble>(data3, &data3->flux_values);
51 flux_divs_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->flux_divs);
52
54 boost::shared_ptr<VectorDouble>(data3, &data3->mass_values);
55
56 mass_dots_ptr3 = boost::shared_ptr<VectorDouble>(data3, &data3->mass_dots);
57 }
58
59 // RDProblem(const int order) : order(order){}
61
62 double global_error;
63
64private:
66 MoFEMErrorCode add_fe(std::string mass_field, std::string flux_field);
67 MoFEMErrorCode set_blockData(std::map<int, BlockData> &block_data_map);
68 MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL, std::string internal);
69 MoFEMErrorCode extract_initial_ents(int block_id, Range &surface);
70 MoFEMErrorCode update_slow_rhs(std::string mass_fiedl,
71 boost::shared_ptr<VectorDouble> &mass_ptr);
72 MoFEMErrorCode push_slow_rhs(std::string mass_field, std::string flux_field,
73 boost::shared_ptr<PreviousData> &data);
74 MoFEMErrorCode update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
75 boost::shared_ptr<PreviousData> &data);
77 update_stiff_rhs(std::string mass_field, std::string flux_field,
78 boost::shared_ptr<VectorDouble> &mass_ptr,
79 boost::shared_ptr<MatrixDouble> &flux_ptr,
80 boost::shared_ptr<VectorDouble> &mass_dot_ptr,
81 boost::shared_ptr<VectorDouble> &flux_div_ptr);
82 MoFEMErrorCode push_stiff_rhs(std::string mass_field, std::string flux_field,
83 boost::shared_ptr<PreviousData> &data,
84 std::map<int, BlockData> &block_map);
85 MoFEMErrorCode update_stiff_lhs(std::string mass_fiedl,
86 std::string flux_field,
87 boost::shared_ptr<VectorDouble> &mass_ptr,
88 boost::shared_ptr<MatrixDouble> &flux_ptr);
89 MoFEMErrorCode push_stiff_lhs(std::string mass_field, std::string flux_field,
90 boost::shared_ptr<PreviousData> &data,
91 std::map<int, BlockData> &block_map);
92
94 MoFEMErrorCode apply_IC(std::string mass_field, Range &surface,
95 boost::shared_ptr<FaceEle> &initial_ele);
96 MoFEMErrorCode apply_BC(std::string flux_field);
98 MoFEMErrorCode post_proc_fields(std::string mass_field,
99 std::string flux_field);
102
107
109 Range natural_bdry_ents;
110
112
113 Range inner_surface1; // nb_species times
114 Range inner_surface2;
115 Range inner_surface3;
116
117 MPI_Comm cOmm;
118 const int rAnk;
119
120 int order;
121 int nb_species;
122
123 std::map<int, BlockData> material_blocks;
124
125 boost::shared_ptr<FaceEle> vol_ele_slow_rhs;
126 boost::shared_ptr<FaceEle> vol_ele_stiff_rhs;
127 boost::shared_ptr<FaceEle> vol_ele_stiff_lhs;
128 boost::shared_ptr<BoundaryEle> natural_bdry_ele_slow_rhs;
129
130 boost::shared_ptr<PostProcFaceOnRefinedMesh> post_proc;
131 boost::shared_ptr<Monitor> monitor_ptr;
132
133 boost::shared_ptr<PreviousData> data1; // nb_species times
134 boost::shared_ptr<PreviousData> data2;
135 boost::shared_ptr<PreviousData> data3;
136
137 boost::shared_ptr<MatrixDouble> flux_values_ptr1; // nb_species times
138 boost::shared_ptr<MatrixDouble> flux_values_ptr2;
139 boost::shared_ptr<MatrixDouble> flux_values_ptr3;
140
141 boost::shared_ptr<VectorDouble> flux_divs_ptr1; // nb_species times
142 boost::shared_ptr<VectorDouble> flux_divs_ptr2;
143 boost::shared_ptr<VectorDouble> flux_divs_ptr3;
144
145 boost::shared_ptr<VectorDouble> mass_values_ptr1; // nb_species times
146 boost::shared_ptr<VectorDouble> mass_values_ptr2;
147 boost::shared_ptr<VectorDouble> mass_values_ptr3;
148
149 boost::shared_ptr<VectorDouble> mass_dots_ptr1; // nb_species times
150 boost::shared_ptr<VectorDouble> mass_dots_ptr2;
151 boost::shared_ptr<VectorDouble> mass_dots_ptr3;
152
153 boost::shared_ptr<ForcesAndSourcesCore> null;
154};
155
156const double ramp_t = 1.0;
157const double sml = 0.0;
158const double T = M_PI / 2.0;
159
161 double operator() (const double x) const {
162 return 1 - abs(x);
163 }
164};
165
167 double operator()(const double x) const {
168 if(x > 0)
169 {return -1;}
170 else
171 {return 1;}
172}
173};
174
176 double operator()(const double x, const double y, const double t) const {
177 double g = cos(T * x) * cos(T * y);
178 if (t <= ramp_t) {
179 return g * t;
180 } else {
181 return g * ramp_t;
182 }
183 }
184};
185
187 FTensor::Tensor1<double, 3> operator()(const double x, const double y,
188 const double t) const {
190 double mx = - T * sin(T * x) * cos(T * y);
191 double my = - T * cos(T * x) * sin(T * y);
192
193
194 if (t <= ramp_t) {
195 grad(0) = mx * t;
196 grad(1) = my * t;
197 } else {
198 grad(0) = mx * ramp_t;
199 grad(1) = my * ramp_t;
200 }
201 grad(2) = 0.0;
202 return grad;
203 }
204};
205
207 double operator()(const double x, const double y, const double t) const {
208 double glap = -2.0 * pow(T, 2) * cos(T * x) * cos(T * y);
209 if (t <= ramp_t) {
210 return glap * t;
211 } else {
212 return glap * ramp_t;
213 }
214 }
215};
216
218 double operator()(const double x, const double y, const double t) const {
219 double gdot = cos(T * x) * cos(T * y);
220 if (t <= ramp_t) {
221 return gdot;
222 } else {
223 return 0;
224 }
225 }
226};
227
234}
235
236MoFEMErrorCode RDProblem::add_fe(std::string mass_field,
237 std::string flux_field) {
241
244
248 1);
249
250 CHKERR simple_interface->setFieldOrder(mass_field, order - 1);
252 CHKERR simple_interface->setFieldOrder("ERROR", 0); // approximation order for error
253
255}
256MoFEMErrorCode RDProblem::set_blockData(std::map<int, BlockData> &block_map) {
259 string name = it->getName();
260 const int id = it->getMeshsetId();
261 if (name.compare(0, 14, "REGION1") == 0) {
262 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
263 id, BLOCKSET, 2, block_map[id].block_ents, true);
264 block_map[id].B0 = 1e-4;
265 block_map[id].block_id = id;
266 } else if (name.compare(0, 14, "REGION2") == 0) {
267 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
268 id, BLOCKSET, 2, block_map[id].block_ents, true);
269 block_map[id].B0 = 1e-4;
270 block_map[id].block_id = id;
271 } else if (name.compare(0, 14, "REGION3") == 0) {
272 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
273 id, BLOCKSET, 2, block_map[id].block_ents, true);
274 block_map[id].B0 = 1e-4;
275 block_map[id].block_id = id;
276 } else if (name.compare(0, 14, "REGION4") == 0) {
277 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
278 id, BLOCKSET, 2, block_map[id].block_ents, true);
279 block_map[id].B0 = 1e-4;
280 block_map[id].block_id = id;
281 } else if (name.compare(0, 14, "REGION5") == 0) {
282 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
283 id, BLOCKSET, 2, block_map[id].block_ents, true);
284 block_map[id].B0 = 1e-4;
285 block_map[id].block_id = id;
286 }
287 }
289}
290
292 std::string natural,
293 std::string internal) {
296 string name = it->getName();
297 if (name.compare(0, 14, natural) == 0) {
298
299 CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
300 natural_bdry_ents, true);
301 } else if (name.compare(0, 14, essential) == 0) {
302 CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
303 essential_bdry_ents, true);
304 } else if (name.compare(0, 14, internal) == 0) {
305 CHKERR it->getMeshsetIdEntitiesByDimension(m_field.get_moab(), 1,
306 internal_edge_ents, true);
307 }
308 }
310}
311
312MoFEMErrorCode RDProblem::extract_initial_ents(int block_id, Range &surface) {
315 BLOCKSET)) {
316 CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
317 block_id, BLOCKSET, 2, surface, true);
318 }
320}
322RDProblem::update_slow_rhs(std::string mass_field,
323 boost::shared_ptr<VectorDouble> &mass_ptr) {
325 vol_ele_slow_rhs->getOpPtrVector().push_back(
326 new OpCalculateScalarFieldValues(mass_field, mass_ptr));
328}
329
331 std::string flux_field,
332 boost::shared_ptr<PreviousData> &data) {
334
335 vol_ele_slow_rhs->getOpPtrVector().push_back(
337
338 natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
339 new OpSkeletonSource(mass_field,
340 KinkFunction(),
341 ExactFunction(),
343
344 natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
346
347 natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
348 new OpEssentialBC(flux_field, essential_bdry_ents));
350}
351
352MoFEMErrorCode RDProblem::update_vol_fe(boost::shared_ptr<FaceEle> &vol_ele,
353 boost::shared_ptr<PreviousData> &data) {
355 auto det_ptr = boost::make_shared<VectorDouble>();
356 auto jac_ptr = boost::make_shared<MatrixDouble>();
357 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
358 vol_ele->getOpPtrVector().push_back(new OpCalculateHOJacForFace(jac_ptr));
359 vol_ele->getOpPtrVector().push_back(
360 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
361 vol_ele->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
362 vol_ele->getOpPtrVector().push_back(
364 vol_ele->getOpPtrVector().push_back(new OpSetInvJacHcurlFace(inv_jac_ptr));
366}
368RDProblem::update_stiff_rhs(std::string mass_field, std::string flux_field,
369 boost::shared_ptr<VectorDouble> &mass_ptr,
370 boost::shared_ptr<MatrixDouble> &flux_ptr,
371 boost::shared_ptr<VectorDouble> &mass_dot_ptr,
372 boost::shared_ptr<VectorDouble> &flux_div_ptr) {
373
375
376 vol_ele_stiff_rhs->getOpPtrVector().push_back(
377 new OpCalculateScalarFieldValues(mass_field, mass_ptr));
378
379 vol_ele_stiff_rhs->getOpPtrVector().push_back(
380 new OpCalculateHVecVectorField<3>(flux_field, flux_ptr));
381
382 vol_ele_stiff_rhs->getOpPtrVector().push_back(
383 new OpCalculateScalarValuesDot(mass_field, mass_dot_ptr));
384
385 vol_ele_stiff_rhs->getOpPtrVector().push_back(
386 new OpCalculateHdivVectorDivergence<3, 2>(flux_field, flux_div_ptr));
388}
389
391 std::string flux_field,
392 boost::shared_ptr<PreviousData> &data,
393 std::map<int, BlockData> &block_map) {
395 vol_ele_stiff_rhs->getOpPtrVector().push_back(
396 new OpAssembleStiffRhsTau<3>(flux_field, data, block_map));
397
398 vol_ele_stiff_rhs->getOpPtrVector().push_back(
402}
403
405RDProblem::update_stiff_lhs(std::string mass_field, std::string flux_field,
406 boost::shared_ptr<VectorDouble> &mass_ptr,
407 boost::shared_ptr<MatrixDouble> &flux_ptr) {
409 vol_ele_stiff_lhs->getOpPtrVector().push_back(
410 new OpCalculateScalarFieldValues(mass_field, mass_ptr));
411
412 vol_ele_stiff_lhs->getOpPtrVector().push_back(
413 new OpCalculateHVecVectorField<3>(flux_field, flux_ptr));
415}
416
418 std::string flux_field,
419 boost::shared_ptr<PreviousData> &data,
420 std::map<int, BlockData> &block_map) {
422 vol_ele_stiff_lhs->getOpPtrVector().push_back(
423 new OpAssembleLhsTauTau<3>(flux_field, data, block_map));
424
425 vol_ele_stiff_lhs->getOpPtrVector().push_back(
426 new OpAssembleLhsVV(mass_field));
427
428 vol_ele_stiff_lhs->getOpPtrVector().push_back(
429 new OpAssembleLhsTauV<3>(flux_field, mass_field, data, block_map));
430
431 vol_ele_stiff_lhs->getOpPtrVector().push_back(
432 new OpAssembleLhsVTau(mass_field, flux_field));
434}
435
438 auto vol_rule = [](int, int, int p) -> int { return 2 * p; };
439 vol_ele_slow_rhs->getRuleHook = vol_rule;
440 natural_bdry_ele_slow_rhs->getRuleHook = vol_rule;
441
442 vol_ele_stiff_rhs->getRuleHook = vol_rule;
443
444 vol_ele_stiff_lhs->getRuleHook = vol_rule;
446}
447
448MoFEMErrorCode RDProblem::apply_IC(std::string mass_field, Range &surface,
449 boost::shared_ptr<FaceEle> &initial_ele) {
451 initial_ele->getOpPtrVector().push_back(
452 new OpInitialMass(mass_field, surface));
454}
455
456MoFEMErrorCode RDProblem::apply_BC(std::string flux_field) {
458 CHKERR m_field.getInterface<ProblemsManager>()->removeDofsOnEntities(
459 "SimpleProblem", flux_field, essential_bdry_ents);
460
462}
465
466 // Vec R, U;
467 // CHKERR DMCreateGlobalVector(dm, &R);
468
469 // CHKERR VecDuplicate(R, &U);
470
471 // DM dm_sub_ff, dm_sub_fm;
472 // ublas::matrix<Mat> nested_matrices(2, 2);
473 // ublas::vector<IS> nested_is(2);
474
475 // CHKERR DMCreate(PETSC_COMM_WORLD, &dm_sub_ff);
476 // CHKERR DMSetType(dm_sub_ff, "DMMOFEM");
477 // CHKERR DMMoFEMCreateSubDM(dm_sub_hh, dm, "SUB_FF");
478
479 // CHKERR DMMoFEMSetSquareProblem(dm_sub_ff, PETSC_TRUE);
480 // CHKERR DMMoFEMAddSubFieldRow(dm_sub_ff, "FLUX1");
481 // CHKERR DMMoFEMAddSubFieldCol(dm_sub_ff, "FLUX1");
482 // CHKERR DMMoFEMAddElement(dm_sub_ff,
483 // simple_interface->getDomainFEName().c_str());
484 // CHKERR DMSetUp(dm_sub_ff);
485
486 // CHKERR DMMoFEMGetSubRowIS(dm_sub_ff, &nested_is[0]);
487 // CHKERR DMCreateMatrix(dm_sub_ff, &nested_matrices(0, 0));
488
489 // vol_ele_stiff_lhs->getFEMethod()->ts_B = nested_matrices(0, 0);
490 // CHKERR TSSetType(ts, TSARKIMEX);
491 // CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
492
493 // CHKERR DMMoFEMTSSetIJacobian(dm_sub_ff, simple_interface->getDomainFEName(),
494 // vol_ele_stiff_lhs, null, null);
495
496 // CHKERR MatAssemblyBegin(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
497 // CHKERR MatAssemblyEnd(nested_matrices(0, 0), MAT_FINAL_ASSEMBLY);
498
499 CHKERR TSSetType(ts, TSARKIMEX);
500 CHKERR TSARKIMEXSetType(ts, TSARKIMEXA2);
501
504
505
508
513
514
516}
517
519 std::string flux_field) {
521 post_proc->addFieldValuesPostProc(mass_field);
522 post_proc->addFieldValuesPostProc(flux_field);
524}
525
531}
534 // Create solution vector
537 CHKERR DMoFEMMeshToLocalVector(dm, X, INSERT_VALUES, SCATTER_FORWARD);
538 // Solve problem
539 double ftime = 1;
540 CHKERR TSSetDM(ts, dm);
541 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
542 CHKERR TSSetSolution(ts, X);
543 CHKERR TSSetFromOptions(ts);
544
545 if (1) {
546 SNES snes;
547 CHKERR TSGetSNES(ts, &snes);
548 KSP ksp;
549 CHKERR SNESGetKSP(snes, &ksp);
550 PC pc;
551 CHKERR KSPGetPC(ksp, &pc);
552 PetscBool is_pcfs = PETSC_FALSE;
553 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
554 // Set up FIELDSPLIT
555 // Only is user set -pc_type fieldsplit
556 if (is_pcfs == PETSC_TRUE) {
557 IS is_mass, is_flux;
558 const MoFEM::Problem *problem_ptr;
559 CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
560 CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
561 problem_ptr->getName(), ROW, "MASS1", 0, 1, &is_mass);
562 CHKERR m_field.getInterface<ISManager>()->isCreateProblemFieldAndRank(
563 problem_ptr->getName(), ROW, "FLUX1", 0, 1, &is_flux);
564 // CHKERR ISView(is_flux, PETSC_VIEWER_STDOUT_SELF);
565 // CHKERR ISView(is_mass, PETSC_VIEWER_STDOUT_SELF);
566
567 CHKERR PCFieldSplitSetIS(pc, NULL, is_flux);
568 CHKERR PCFieldSplitSetIS(pc, NULL, is_mass);
569
570
571 CHKERR ISDestroy(&is_flux);
572 CHKERR ISDestroy(&is_mass);
573 }
574 }
575
576 CHKERR TSSolve(ts, X);
578}
579
582 global_error = 0;
583 // set nb_species
584 CHKERR setup_system(); // only once
585 nb_species = nb_sp;
586 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
587 CHKERR add_fe("MASS1", "FLUX1"); // nb_species times
588 if (nb_species == 2 || nb_species == 3) {
589 CHKERR add_fe("MASS2", "FLUX2");
590 if (nb_species == 3) {
591 CHKERR add_fe("MASS3", "FLUX3");
592 }
593 }
594 }
595
597
599
600 CHKERR extract_bd_ents("ESSENTIAL", "NATURAL", "INTERNAL"); // nb_species times
601
602 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
605 if (nb_species == 1) {
606 vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
607 "MASS1", data1, data1, data1, material_blocks));
608 } else if (nb_species == 2 || nb_species == 3) {
611 if (nb_species == 2) {
612 vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
613 "MASS1", data1, data2, data2, material_blocks));
614 } else if (nb_species == 3) {
617 vol_ele_slow_rhs->getOpPtrVector().push_back(new OpComputeSlowValue(
618 "MASS1", data1, data2, data3, material_blocks));
619 }
620 }
621 }
622 natural_bdry_ele_slow_rhs->getOpPtrVector().push_back(
624
625 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
626 CHKERR push_slow_rhs("MASS1", "FLUX1", data1); // nb_species times
627 if (nb_species == 2 || nb_species == 3) {
628 CHKERR push_slow_rhs("MASS2", "FLUX2", data2);
629 if (nb_species == 3) {
630 CHKERR push_slow_rhs("MASS3", "FLUX3", data3);
631 }
632 }
633 }
634
636
637 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
638 CHKERR update_stiff_rhs("MASS1", "FLUX1", mass_values_ptr1,
640 CHKERR push_stiff_rhs("MASS1", "FLUX1", data1,
641 material_blocks); // nb_species times
642 if (nb_species == 2 || nb_species == 3) {
643 CHKERR update_stiff_rhs("MASS2", "FLUX2", mass_values_ptr2,
645 CHKERR push_stiff_rhs("MASS2", "FLUX2", data2, material_blocks);
646 if (nb_species == 3) {
647 CHKERR update_stiff_rhs("MASS3", "FLUX3", mass_values_ptr3,
650 CHKERR push_stiff_rhs("MASS3", "FLUX3", data3, material_blocks);
651 }
652 }
653 }
654
656
657 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
658 CHKERR update_stiff_lhs("MASS1", "FLUX1", mass_values_ptr1,
660 CHKERR push_stiff_lhs("MASS1", "FLUX1", data1,
661 material_blocks); // nb_species times
662 if (nb_species == 2 || nb_species == 3) {
663 CHKERR update_stiff_lhs("MASS2", "FLUX2", mass_values_ptr2,
665 CHKERR push_stiff_lhs("MASS2", "FLUX2", data2, material_blocks);
666 if (nb_species == 3) {
667 CHKERR update_stiff_lhs("MASS3", "FLUX3", mass_values_ptr3,
669 CHKERR push_stiff_lhs("MASS3", "FLUX3", data3, material_blocks);
670 }
671 }
672 }
673
674
675 // vol_ele_stiff_lhs->getOpPtrVector().push_back(
676 // new OpError(ExactFunction(), ExactFunctionLap(), ExactFunctionGrad(), data1, material_blocks, global_error));
680 boost::shared_ptr<FaceEle> initial_mass_ele(new FaceEle(m_field));
681
682 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
684 initial_mass_ele); // nb_species times
685 if (nb_species == 2 || nb_species == 3) {
686 CHKERR apply_IC("MASS2", inner_surface2, initial_mass_ele);
687 if (nb_species == 3) {
688 CHKERR apply_IC("MASS3", inner_surface3, initial_mass_ele);
689 }
690 }
691 }
693 initial_mass_ele);
694
695 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
696 CHKERR apply_BC("FLUX1"); // nb_species times
697 if (nb_species == 2 || nb_species == 3) {
698 CHKERR apply_BC("FLUX2");
699 if (nb_species == 3) {
700 CHKERR apply_BC("FLUX3");
701 }
702 }
703 }
704
705 CHKERR loop_fe(); // only once
706 post_proc->generateReferenceElementMesh(); // only once
707
708
709
710 if (nb_species == 1 || nb_species == 2 || nb_species == 3) {
711 CHKERR post_proc_fields("MASS1", "FLUX1");
712 post_proc->addFieldValuesPostProc("ERROR");
713 if (nb_species == 2 || nb_species == 3) {
714 CHKERR post_proc_fields("MASS2", "FLUX2");
715 if (nb_species == 3) {
716 CHKERR post_proc_fields("MASS3", "FLUX3");
717 }
718 }
719 }
720
721 //the double global_vector is a global_vector which sums errors in each element
722
723 // Vec gVec;
724 // CHKERR VecCreateMPI(PETSC_COMM_WORLD, 1, 1, &gVec);
725 // double err_per_step = 0;
726
727 // auto compute_global_error = [&gVec](double &g_err, double &err_out) {
728 // MoFEMFunctionBegin;
729
730 // CHKERR VecZeroEntries(gVec);
731 // CHKERR VecAssemblyBegin(gVec);
732 // CHKERR VecAssemblyEnd(gVec);
733
734 // int ind[1] = {0};
735 // CHKERR VecSetValues(gVec, 1, ind, &err_out, ADD_VALUES);
736 // CHKERR VecAssemblyBegin(gVec);
737 // CHKERR VecAssemblyEnd(gVec);
738
739 // CHKERR VecGetValues(gVec, 1, ind, &err_out);
740
741 // MoFEMFunctionReturn(0);
742 // };
743
744 // Vec error_per_proc;
745 // CHKERR VecCreateMPI(cOmm, 1, PETSC_DECIDE, &error_per_proc);
746 // auto get_global_error = [&]() {
747 // MoFEMFunctionBegin;
748 // CHKERR VecSetValue(error_per_proc, m_field.get_comm_rank(), global_error, INSERT_VALUES);
749 // MoFEMFunctionReturn(0);
750 // };
751
752 // auto reinitialize_global = [&]() {
753 // MoFEMFunctionBegin;
754 // CHKERR get_global_error();
755 // CHKERR VecAssemblyBegin(error_per_proc);
756 // CHKERR VecAssemblyEnd(error_per_proc);
757 // double error_sum;
758 // CHKERR VecSum(error_per_proc, &error_sum);
759 // CHKERR PetscPrintf(PETSC_COMM_WORLD, "Error : %3.4e \n",
760 // error_sum);
761 // global_error = 0;
762 // MoFEMFunctionReturn(0);
763 // };
764
765 // vol_ele_stiff_lhs->postProcessHook = reinitialize_global;
766
767
768
769 monitor_ptr = boost::shared_ptr<Monitor>(
770 new Monitor(cOmm, rAnk, dm, post_proc, global_error)); // nb_species times
771 CHKERR output_result(); // only once
772 CHKERR solve(); // only once
774}
775
776int main(int argc, char *argv[]) {
777 const char param_file[] = "param_file.petsc";
779 try {
780 moab::Core mb_instance;
781 moab::Interface &moab = mb_instance;
782 MoFEM::Core core(moab);
783 DMType dm_name = "DMMOFEM";
784 CHKERR DMRegister_MoFEM(dm_name);
785
786 int order = 1;
787 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
788 int nb_species = 1;
789 RDProblem reac_diff_problem(core, order+1);
790 CHKERR reac_diff_problem.run_analysis(nb_species);
791 }
794 return 0;
795}
static Index< 'p', 3 > p
std::string param_file
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:79
@ L2
field with C-1 continuity
Definition: definitions.h:101
@ HCURL
field with continuous tangents
Definition: definitions.h:99
#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 DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:385
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.
int main(int argc, char *argv[])
static char help[]
const double sml
const double ramp_t
const double T
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
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
OpSetContravariantPiolaTransformOnEdge2D OpSetContrariantPiolaTransformOnEdge
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
constexpr double t
plate stiffness
Definition: plate.cpp:76
constexpr double g
double operator()(const double x) const
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
double operator()(const double x) 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.
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:33
Interface for managing meshsets containing materials and boundary conditions.
Get vector field for H-div approximation.
Calculate divergence of vector field.
Get value at integration points for scalar field.
Make Hdiv space from Hcurl space in 2d.
keeps basic data about problem
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:35
const std::string getBoundaryFEName() const
Get the Boundary FE Name.
Definition: Simple.hpp:328
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 addBoundaryField(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 boundary.
Definition: Simple.cpp:396
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.
MoFEMErrorCode update_slow_rhs(std::string mass_fiedl, boost::shared_ptr< VectorDouble > &mass_ptr)
boost::shared_ptr< MatrixDouble > flux_values_ptr2
boost::shared_ptr< FaceEle > vol_ele_stiff_lhs
boost::shared_ptr< VectorDouble > flux_divs_ptr3
std::map< int, BlockData > material_blocks
SmartPetscObj< TS > ts
MoFEMErrorCode update_stiff_rhs()
Range inner_surface1
SmartPetscObj< DM > dm
MoFEMErrorCode apply_BC(std::string flux_field)
boost::shared_ptr< Monitor > monitor_ptr
MoFEMErrorCode add_fe()
MoFEMErrorCode push_slow_rhs()
boost::shared_ptr< PreviousData > data1
boost::shared_ptr< VectorDouble > mass_values_ptr3
MoFEMErrorCode update_vol_fe(boost::shared_ptr< FaceEle > &vol_ele, boost::shared_ptr< PreviousData > &data)
Range inner_surface2
Range essential_bdry_ents
MoFEMErrorCode post_proc_fields()
boost::shared_ptr< ForcesAndSourcesCore > null
boost::shared_ptr< VectorDouble > flux_divs_ptr2
MoFEMErrorCode loop_fe()
MoFEM::Interface & m_field
Simple * simple_interface
Range inner_surface3
MoFEMErrorCode run_analysis()
MPI_Comm cOmm
boost::shared_ptr< MatrixDouble > flux_values_ptr3
boost::shared_ptr< VectorDouble > flux_divs_ptr1
MoFEMErrorCode extract_bd_ents(std::string ESSENTIAL, std::string NATURAL)
Range internal_edge_ents
double global_error
boost::shared_ptr< BoundaryEle > natural_bdry_ele_slow_rhs
MoFEMErrorCode set_blockData(std::map< int, BlockData > &block_data_map)
const int rAnk
RDProblem(MoFEM::Core &core, const int order)
MoFEMErrorCode push_stiff_lhs()
MoFEMErrorCode push_stiff_rhs()
MoFEMErrorCode extract_initial_ents(int block_id, Range &surface)
MoFEMErrorCode solve()
boost::shared_ptr< VectorDouble > mass_dots_ptr3
boost::shared_ptr< VectorDouble > mass_values_ptr2
MoFEMErrorCode apply_IC(std::string mass_field, Range &surface, boost::shared_ptr< FaceEle > &initial_ele, double &init_val)
Range natural_bdry_ents
boost::shared_ptr< VectorDouble > mass_dots_ptr1
boost::shared_ptr< VectorDouble > mass_values_ptr1
MoFEMErrorCode update_stiff_lhs(std::string mass_fiedl, std::string flux_field, boost::shared_ptr< VectorDouble > &mass_ptr, boost::shared_ptr< MatrixDouble > &flux_ptr)
boost::shared_ptr< VectorDouble > mass_dots_ptr2
boost::shared_ptr< PreviousData > data2
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()
boost::shared_ptr< PreviousData > data3
MoFEMErrorCode update_stiff_rhs(std::string mass_field, std::string flux_field, boost::shared_ptr< VectorDouble > &mass_ptr, boost::shared_ptr< MatrixDouble > &flux_ptr, boost::shared_ptr< VectorDouble > &mass_dot_ptr, boost::shared_ptr< VectorDouble > &flux_div_ptr)
MoFEMErrorCode set_integration_rule()
boost::shared_ptr< PostProcFaceOnRefinedMesh > post_proc
boost::shared_ptr< MatrixDouble > flux_values_ptr1
boost::shared_ptr< FaceEle > vol_ele_stiff_rhs