v0.13.1
ElecPhysOperators2D.hpp
Go to the documentation of this file.
1#ifndef __ELECPHYSOPERATORS2D_HPP__
2#define __ELECPHYSOPERATORS2D_HPP__
3
4#include <stdlib.h>
5#include <BasicFiniteElements.hpp>
6
7namespace ElectroPhysiology {
8
9using FaceEle = MoFEM::FaceElementForcesAndSourcesCoreSwitch<0>;
10
11using BoundaryEle = MoFEM::EdgeElementForcesAndSourcesCoreSwitch<0>;
12
15
17
19
20double factor = 1.0/ 12.9;
21
22const double B = 0.0;
23const double B_epsilon = 0.0;
24
25const int save_every_nth_step = 16;
26
27const double essen_value = 0;
28
30
31// problem parameters
32const double alpha = 0.01;
33const double gma = 0.002;
34const double b = 0.15;
35const double c = 8.00;
36const double mu1 = 0.20;
37const double mu2 = 0.30;
38
39struct BlockData {
41
43
44 double B0; // species mobility
45
47 :
48 B0(0.2) {}
49};
50
52
53const double D_tilde = 1e-2;
54
55struct PreviousData {
58
61
62
65
67 jac.resize(2, 2, false);
68 inv_jac.resize(2, 2, false);
69 }
70};
71struct OpEssentialBC : public OpBoundaryEle {
72 OpEssentialBC(const std::string &flux_field, Range &essential_bd_ents)
73 : OpBoundaryEle(flux_field, OpBoundaryEle::OPROW),
75
78 int nb_dofs = data.getIndices().size();
79 if (nb_dofs) {
80 EntityHandle fe_ent = getFEEntityHandle();
81 bool is_essential =
82 (essential_bd_ents.find(fe_ent) != essential_bd_ents.end());
83 if (is_essential) {
84 int nb_gauss_pts = getGaussPts().size2();
85 int size2 = data.getN().size2();
86 if (3 * nb_dofs != static_cast<int>(data.getN().size2()))
87 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
88 "wrong number of dofs");
89 nN.resize(nb_dofs, nb_dofs, false);
90 nF.resize(nb_dofs, false);
91 nN.clear();
92 nF.clear();
93
94 auto t_row_tau = data.getFTensor1N<3>();
95
96 auto dir = getDirection();
97 double len = sqrt(dir[0] * dir[0] + dir[1] * dir[1] + dir[2] * dir[2]);
98
99 FTensor::Tensor1<double, 3> t_normal(-dir[1] / len, dir[0] / len,
100 dir[2] / len);
101
102 auto t_w = getFTensor0IntegrationWeight();
103 const double vol = getMeasure();
104
105 for (int gg = 0; gg < nb_gauss_pts; gg++) {
106 const double a = t_w * vol;
107 for (int rr = 0; rr != nb_dofs; rr++) {
108 auto t_col_tau = data.getFTensor1N<3>(gg, 0);
109 nF[rr] += a * essen_value * t_row_tau(i) * t_normal(i);
110 for (int cc = 0; cc != nb_dofs; cc++) {
111 nN(rr, cc) += a * (t_row_tau(i) * t_normal(i)) *
112 (t_col_tau(i) * t_normal(i));
113 ++t_col_tau;
114 }
115 ++t_row_tau;
116 }
117 ++t_w;
118 }
119
121 cholesky_solve(nN, nF, ublas::lower());
122
123 for (auto &dof : data.getFieldDofs()) {
124 dof->getFieldData() = nF[dof->getEntDofIdx()];
125 }
126 }
127 }
129 }
130
131private:
134 Range &essential_bd_ents;
135};
136
137// Assembly of system mass matrix
138// //***********************************************
139
140// Mass matrix corresponding to the flux equation.
141// 01. Note that it is an identity matrix
142struct OpInitialMass : public OpFaceEle {
143 OpInitialMass(const std::string &mass_field, Range &inner_surface, double &init_val)
144 : OpFaceEle(mass_field, OpFaceEle::OPROW), innerSurface(inner_surface), initVal(init_val) {}
147 double &initVal;
148 Range &innerSurface;
151 int nb_dofs = data.getFieldData().size();
152 if (nb_dofs) {
153 EntityHandle fe_ent = getFEEntityHandle();
154 bool is_inner_side = (innerSurface.find(fe_ent) != innerSurface.end());
155 if (is_inner_side) {
156 int nb_gauss_pts = getGaussPts().size2();
157 if (nb_dofs != static_cast<int>(data.getN().size2()))
158 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
159 "wrong number of dofs");
160 nN.resize(nb_dofs, nb_dofs, false);
161 nF.resize(nb_dofs, false);
162 nN.clear();
163 nF.clear();
164
165 auto t_row_mass = data.getFTensor0N();
166 auto t_w = getFTensor0IntegrationWeight();
167 const double vol = getMeasure();
168
169 for (int gg = 0; gg < nb_gauss_pts; gg++) {
170 const double a = t_w * vol;
171 double r = initVal;
172 for (int rr = 0; rr != nb_dofs; rr++) {
173 auto t_col_mass = data.getFTensor0N(gg, 0);
174 nF[rr] += a * r * t_row_mass;
175 for (int cc = 0; cc != nb_dofs; cc++) {
176 nN(rr, cc) += a * t_row_mass * t_col_mass;
177 ++t_col_mass;
178 }
179 ++t_row_mass;
180 }
181 ++t_w;
182 }
183
185 cholesky_solve(nN, nF, ublas::lower());
186
187 for (auto &dof : data.getFieldDofs()) {
188 dof->getFieldData() = nF[dof->getEntDofIdx()];
189
190 // this is only to check
191 // data.getFieldData()[dof->getEntDofIdx()] = nF[dof->getEntDofIdx()];
192 }
193 }
194 }
196 }
197};
198
199struct OpSolveRecovery : public OpFaceEle {
200 typedef boost::function<double(const double, const double, const double)>
202 OpSolveRecovery(const std::string &mass_field,
203 boost::shared_ptr<PreviousData> &data_u,
204 boost::shared_ptr<PreviousData> &data_v, Method runge_kutta4)
205 : OpFaceEle(mass_field, OpFaceEle::OPROW), dataU(data_u), dataV(data_v),
206 rungeKutta4(runge_kutta4) {}
207 boost::shared_ptr<PreviousData> dataU;
208 boost::shared_ptr<PreviousData> dataV;
210
213 double initVal;
216 int nb_dofs = data.getFieldData().size();
217 if (nb_dofs) {
218 int nb_gauss_pts = getGaussPts().size2();
219 if (nb_dofs != static_cast<int>(data.getN().size2()))
220 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
221 "wrong number of dofs");
222 nN.resize(nb_dofs, nb_dofs, false);
223 nF.resize(nb_dofs, false);
224 nN.clear();
225 nF.clear();
226
227 auto t_val_u = getFTensor0FromVec(dataU->mass_values);
228 auto t_val_v = getFTensor0FromVec(dataV->mass_values);
229
230 double dt;
231 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
232
233 auto t_row_mass = data.getFTensor0N();
234 auto t_w = getFTensor0IntegrationWeight();
235 const double vol = getMeasure();
236
237 for (int gg = 0; gg < nb_gauss_pts; gg++) {
238 const double a = t_w * vol;
239 const double vn = rungeKutta4(t_val_u, t_val_v, dt);
240 for (int rr = 0; rr != nb_dofs; rr++) {
241 auto t_col_mass = data.getFTensor0N(gg, 0);
242 nF[rr] += a * vn * t_row_mass;
243 for (int cc = 0; cc != nb_dofs; cc++) {
244 nN(rr, cc) += a * t_row_mass * t_col_mass;
245 ++t_col_mass;
246 }
247 ++t_row_mass;
248 }
249 ++t_w;
250 ++t_val_u;
251 ++t_val_v;
252 }
253
255 cholesky_solve(nN, nF, ublas::lower());
256
257 for (auto &dof : data.getFieldDofs()) {
258 dof->getFieldData() = nF[dof->getEntDofIdx()];
259
260 // this is only to check
261 // data.getFieldData()[dof->getEntDofIdx()] = nF[dof->getEntDofIdx()];
262 }
263 }
265 }
266};
267
268// Assembly of RHS for explicit (slow)
269// part//**************************************
270
271// 2. RHS for explicit part of the mass balance equation
272struct OpAssembleSlowRhsV : OpFaceEle // R_V
273{
274 typedef boost::function<double(const double, const double)>
276 OpAssembleSlowRhsV(std::string mass_field,
277 boost::shared_ptr<PreviousData> &common_datau,
278 boost::shared_ptr<PreviousData> &common_datav, FUVal rhs_u)
279 : OpFaceEle(mass_field, OpFaceEle::OPROW), commonDatau(common_datau),
280 commonDatav(common_datav), rhsU(rhs_u) {}
281
284 // cerr << "In OpAssembleSlowRhsV...." << endl;
285 const int nb_dofs = data.getIndices().size();
286 if (nb_dofs) {
287 // cerr << "In SlowRhsV..." << endl;
288 if (nb_dofs != static_cast<int>(data.getN().size2()))
289 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
290 "wrong number of dofs");
291 vecF.resize(nb_dofs, false);
292 mat.resize(nb_dofs, nb_dofs, false);
293 vecF.clear();
294 mat.clear();
295 const int nb_integration_pts = getGaussPts().size2();
296 auto t_u_value = getFTensor0FromVec(commonDatau->mass_values);
297 auto t_v_value = getFTensor0FromVec(commonDatav->mass_values);
298 auto t_row_v_base = data.getFTensor0N();
299 auto t_w = getFTensor0IntegrationWeight();
300 const double vol = getMeasure();
301
302 const double ct = getFEMethod()->ts_t - 0.01;
303 auto t_coords = getFTensor1CoordsAtGaussPts();
304 for (int gg = 0; gg != nb_integration_pts; ++gg) {
305 const double a = vol * t_w;
306
307 // double u_dot = exactDot(t_coords(NX), t_coords(NY), ct);
308 // double u_lap = exactLap(t_coords(NX), t_coords(NY), ct);
309
310 // double f = u_dot - u_lap;
311 double f = t_u_value * (1.0 - t_u_value);
312 for (int rr = 0; rr != nb_dofs; ++rr) {
313 double rhs = rhsU(t_u_value, t_v_value);
314 auto t_col_v_base = data.getFTensor0N(gg, 0);
315 vecF[rr] += a * rhs * t_row_v_base;
316 // vecF[rr] += a * f * t_row_v_base;
317 for (int cc = 0; cc != nb_dofs; ++cc) {
318 mat(rr, cc) += a * t_row_v_base * t_col_v_base;
319 ++t_col_v_base;
320 }
321 ++t_row_v_base;
322 }
323 ++t_u_value;
324 ++t_v_value;
325 ++t_w;
326 ++t_coords;
327 }
329 cholesky_solve(mat, vecF, ublas::lower());
330
331 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
332 PETSC_TRUE);
333 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
334 ADD_VALUES);
335 }
337 }
338
339private:
340 boost::shared_ptr<PreviousData> commonDatau;
341 boost::shared_ptr<PreviousData> commonDatav;
345
346
350};
351
352// // 5. RHS contribution of the natural boundary condition
353// struct OpAssembleNaturalBCRhsTau : OpFaceEle // R_tau_2
354// {
355// OpAssembleNaturalBCRhsTau(std::string flux_field, Range &natural_bd_ents)
356// : OpFaceEle(flux_field, OpFaceEle::OPROW),
357// natural_bd_ents(natural_bd_ents) {}
358
359// MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
360// MoFEMFunctionBegin;
361// const int nb_dofs = data.getIndices().size();
362
363// if (nb_dofs) {
364// EntityHandle row_side_ent = data.getFieldDofs()[0]->getEnt();
365
366// bool is_natural =
367// (natural_bd_ents.find(row_side_ent) != natural_bd_ents.end());
368// if (is_natural) {
369// // cerr << "In NaturalBCRhsTau..." << endl;
370// vecF.resize(nb_dofs, false);
371// vecF.clear();
372// const int nb_integration_pts = getGaussPts().size2();
373// auto t_tau_base = data.getFTensor1N<3>();
374
375// auto dir = getDirection();
376// FTensor::Tensor1<double, 3> t_normal(-dir[1], dir[0], dir[2]);
377
378// auto t_w = getFTensor0IntegrationWeight();
379
380// for (int gg = 0; gg != nb_integration_pts; ++gg) {
381// const double a = t_w;
382// for (int rr = 0; rr != nb_dofs; ++rr) {
383// vecF[rr] += (t_tau_base(i) * t_normal(i) * natu_value) * a;
384// ++t_tau_base;
385// }
386// ++t_w;
387// }
388// CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
389// PETSC_TRUE);
390// CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
391// ADD_VALUES);
392// }
393// }
394// MoFEMFunctionReturn(0);
395// }
396
397// private:
398// VectorDouble vecF;
399// Range natural_bd_ents;
400// };
401
402// Assembly of RHS for the implicit (stiff) part excluding the essential
403// boundary //**********************************
404// 3. Assembly of F_tau excluding the essential boundary condition
405template <int dim>
406struct OpAssembleStiffRhsTau : OpFaceEle // F_tau_1
407{
408 OpAssembleStiffRhsTau(std::string flux_field,
409 boost::shared_ptr<PreviousData> &data,
410 std::map<int, BlockData> &block_map)
411 : OpFaceEle(flux_field, OpFaceEle::OPROW), commonData(data),
413
416
417 const int nb_dofs = data.getIndices().size();
418 if (nb_dofs) {
419 // auto find_block_data = [&]() {
420 // EntityHandle fe_ent = getFEEntityHandle();
421 // BlockData *block_raw_ptr = nullptr;
422 // for (auto &m : setOfBlock) {
423 // if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
424 // block_raw_ptr = &m.second;
425 // break;
426 // }
427 // }
428 // return block_raw_ptr;
429 // };
430
431 // auto block_data_ptr = find_block_data();
432 // if (!block_data_ptr)
433 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
434 // auto &block_data = *block_data_ptr;
435
436 vecF.resize(nb_dofs, false);
437 vecF.clear();
438
439 const int nb_integration_pts = getGaussPts().size2();
440 auto t_flux_value = getFTensor1FromMat<3>(commonData->flux_values);
441 auto t_mass_value = getFTensor0FromVec(commonData->mass_values);
442 auto t_tau_base = data.getFTensor1N<3>();
443
444 auto t_tau_grad = data.getFTensor2DiffN<3, 2>();
445
446 auto t_w = getFTensor0IntegrationWeight();
447 const double vol = getMeasure();
448
449 for (int gg = 0; gg < nb_integration_pts; ++gg) {
450
451 const double K = B_epsilon + (block.B0 + B * t_mass_value);
452 const double K_inv = 1. / K;
453 const double a = vol * t_w;
454 for (int rr = 0; rr < nb_dofs; ++rr) {
455 double div_base = t_tau_grad(0, 0) + t_tau_grad(1, 1);
456 vecF[rr] += (K_inv * t_tau_base(i) * t_flux_value(i) -
457 div_base * t_mass_value) *
458 a;
459 ++t_tau_base;
460 ++t_tau_grad;
461 }
462 ++t_flux_value;
463 ++t_mass_value;
464 ++t_w;
465 }
466 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
467 PETSC_TRUE);
468 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
469 ADD_VALUES);
470 }
472 }
473
474private:
475 boost::shared_ptr<PreviousData> commonData;
477 std::map<int, BlockData> setOfBlock;
478};
479// 4. Assembly of F_v
480template <int dim>
481struct OpAssembleStiffRhsV : OpFaceEle // F_V
482{
483 typedef boost::function<double(const double, const double)>
485 OpAssembleStiffRhsV(std::string flux_field,
486 boost::shared_ptr<PreviousData> &datau,
487 boost::shared_ptr<PreviousData> &datav, FUval rhs_u,
488 std::map<int, BlockData> &block_map, Range &stim_region)
489 : OpFaceEle(flux_field, OpFaceEle::OPROW), commonDatau(datau),
490 commonDatav(datav), setOfBlock(block_map), rhsU(rhs_u),
491 stimRegion(stim_region) {}
492
493 Range &stimRegion;
497 const int nb_dofs = data.getIndices().size();
498 // cerr << "In StiffRhsV ..." << endl;
499 if (nb_dofs) {
500 // auto find_block_data = [&]() {
501 // EntityHandle fe_ent = getFEEntityHandle();
502 // BlockData *block_raw_ptr = nullptr;
503 // for (auto &m : setOfBlock) {
504 // if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
505 // block_raw_ptr = &m.second;
506 // break;
507 // }
508 // }
509 // return block_raw_ptr;
510 // };
511
512 // auto block_data_ptr = find_block_data();
513 // if (!block_data_ptr)
514 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
515 // auto &block_data = *block_data_ptr;
516
517 vecF.resize(nb_dofs, false);
518 vecF.clear();
519 const int nb_integration_pts = getGaussPts().size2();
520 auto t_u_value = getFTensor0FromVec(commonDatau->mass_values);
521 auto t_v_value = getFTensor0FromVec(commonDatav->mass_values);
522
523 auto t_mass_dot = getFTensor0FromVec(commonDatau->mass_dots);
524 auto t_flux_div = getFTensor0FromVec(commonDatau->flux_divs);
525 auto t_row_v_base = data.getFTensor0N();
526 auto t_w = getFTensor0IntegrationWeight();
527 const double vol = getMeasure();
528 auto t_coords = getFTensor1CoordsAtGaussPts();
529
530 const double c_time = getFEMethod()->ts_t;
531
532 double dt;
533 CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
534
535 double stim = 0.0;
536
537 double T = 627.0;
538 double duration = 5.0;
539
540 if (T - dt < c_time && c_time <= T + duration) {
541 EntityHandle stim_ent = getFEEntityHandle();
542 if (stimRegion.find(stim_ent) != stimRegion.end()) {
543 stim = 40.0;
544 } else {
545 stim = 0.0;
546 }
547 }
548 for (int gg = 0; gg < nb_integration_pts; ++gg) {
549 const double a = vol * t_w;
550 // double u_dot = exactDot(t_coords(NX), t_coords(NY), ct);
551 // double u_lap = exactLap(t_coords(NX), t_coords(NY), ct);
552
553 // double f = u_dot - block_data.B0 * u_lap;
554 double rhsu = rhsU(t_u_value, t_v_value);
555 for (int rr = 0; rr < nb_dofs; ++rr) {
556 vecF[rr] += (t_row_v_base * (t_mass_dot + t_flux_div - 0*rhsu - factor * stim)) * a;
557 ++t_row_v_base;
558 }
559 ++t_mass_dot;
560 ++t_flux_div;
561 ++t_u_value;
562 ++t_v_value;
563 ++t_w;
564 ++t_coords;
565 }
566 CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
567 PETSC_TRUE);
568 CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
569 ADD_VALUES);
570 }
572 }
573
574private:
575 boost::shared_ptr<PreviousData> commonDatau;
576 boost::shared_ptr<PreviousData> commonDatav;
578 std::map<int, BlockData> setOfBlock;
579
580 // FVal exactValue;
581 // FVal exactDot;
582 // FVal exactLap;
583
586};
587
588// Tangent operator
589// //**********************************************
590// 7. Tangent assembly for F_tautau excluding the essential boundary condition
591template <int dim>
592struct OpAssembleLhsTauTau : OpFaceEle // A_TauTau_1
593{
594 OpAssembleLhsTauTau(std::string flux_field,
595 boost::shared_ptr<PreviousData> &commonData,
596 std::map<int, BlockData> &block_map)
597 : OpFaceEle(flux_field, flux_field, OpFaceEle::OPROWCOL),
599 sYmm = true;
600 }
601
602 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
603 EntityType col_type, EntData &row_data,
604 EntData &col_data) {
606 const int nb_row_dofs = row_data.getIndices().size();
607 const int nb_col_dofs = col_data.getIndices().size();
608
609 if (nb_row_dofs && nb_col_dofs) {
610 // auto find_block_data = [&]() {
611 // EntityHandle fe_ent = getFEEntityHandle();
612 // BlockData *block_raw_ptr = nullptr;
613 // for (auto &m : setOfBlock) {
614 // if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
615 // block_raw_ptr = &m.second;
616 // break;
617 // }
618 // }
619 // return block_raw_ptr;
620 // };
621
622 // auto block_data_ptr = find_block_data();
623 // if (!block_data_ptr)
624 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
625 // auto &block_data = *block_data_ptr;
626
627 mat.resize(nb_row_dofs, nb_col_dofs, false);
628 mat.clear();
629 const int nb_integration_pts = getGaussPts().size2();
630 auto t_mass_value = getFTensor0FromVec(commonData->mass_values);
631
632 auto t_row_tau_base = row_data.getFTensor1N<3>();
633
634 auto t_w = getFTensor0IntegrationWeight();
635 const double vol = getMeasure();
636
637 for (int gg = 0; gg != nb_integration_pts; ++gg) {
638 const double a = vol * t_w;
639 const double K = B_epsilon + (block.B0 + B * t_mass_value);
640 const double K_inv = 1. / K;
641 for (int rr = 0; rr != nb_row_dofs; ++rr) {
642 auto t_col_tau_base = col_data.getFTensor1N<3>(gg, 0);
643 for (int cc = 0; cc != nb_col_dofs; ++cc) {
644 mat(rr, cc) += (K_inv * t_row_tau_base(i) * t_col_tau_base(i)) * a;
645 ++t_col_tau_base;
646 }
647 ++t_row_tau_base;
648 }
649 ++t_mass_value;
650 ++t_w;
651 }
652 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
653 ADD_VALUES);
654 if (row_side != col_side || row_type != col_type) {
655 transMat.resize(nb_col_dofs, nb_row_dofs, false);
656 noalias(transMat) = trans(mat);
657 CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
658 &transMat(0, 0), ADD_VALUES);
659 }
660 }
662 }
663
664private:
665 boost::shared_ptr<PreviousData> commonData;
668 std::map<int, BlockData> setOfBlock;
669};
670
671// 9. Assembly of tangent for F_tau_v excluding the essential bc
672template <int dim>
673struct OpAssembleLhsTauV : OpFaceEle // E_TauV
674{
675 OpAssembleLhsTauV(std::string flux_field, std::string mass_field,
676 boost::shared_ptr<PreviousData> &data,
677 std::map<int, BlockData> &block_map)
678 : OpFaceEle(flux_field, mass_field, OpFaceEle::OPROWCOL),
680 sYmm = false;
681 }
682
683 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
684 EntityType col_type, EntData &row_data,
685 EntData &col_data) {
687 const int nb_row_dofs = row_data.getIndices().size();
688 const int nb_col_dofs = col_data.getIndices().size();
689
690 if (nb_row_dofs && nb_col_dofs) {
691 // auto find_block_data = [&]() {
692 // EntityHandle fe_ent = getFEEntityHandle();
693 // BlockData *block_raw_ptr = nullptr;
694 // for (auto &m : setOfBlock) {
695 // if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
696 // block_raw_ptr = &m.second;
697 // break;
698 // }
699 // }
700 // return block_raw_ptr;
701 // };
702
703 // auto block_data_ptr = find_block_data();
704 // if (!block_data_ptr)
705 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
706 // auto &block_data = *block_data_ptr;
707 mat.resize(nb_row_dofs, nb_col_dofs, false);
708 mat.clear();
709 const int nb_integration_pts = getGaussPts().size2();
710 auto t_w = getFTensor0IntegrationWeight();
711 auto t_row_tau_base = row_data.getFTensor1N<3>();
712
713 auto t_row_tau_grad = row_data.getFTensor2DiffN<3, 2>();
714 auto t_mass_value = getFTensor0FromVec(commonData->mass_values);
715 auto t_flux_value = getFTensor1FromMat<3>(commonData->flux_values);
716 const double vol = getMeasure();
717
718 for (int gg = 0; gg != nb_integration_pts; ++gg) {
719 const double a = vol * t_w;
720 const double K = B_epsilon + (block.B0 + B * t_mass_value);
721 const double K_inv = 1. / K;
722 const double K_diff = B;
723
724 for (int rr = 0; rr != nb_row_dofs; ++rr) {
725 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
726 for (int cc = 0; cc != nb_col_dofs; ++cc) {
727 double div_row_base = t_row_tau_grad(0, 0) + t_row_tau_grad(1, 1);
728 mat(rr, cc) += (-(t_row_tau_base(i) * t_flux_value(i) * K_inv *
729 K_inv * K_diff * t_col_v_base) -
730 (div_row_base * t_col_v_base)) *
731 a;
732 ++t_col_v_base;
733 }
734 ++t_row_tau_base;
735 ++t_row_tau_grad;
736 }
737 ++t_w;
738 ++t_mass_value;
739 ++t_flux_value;
740 }
741 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
742 ADD_VALUES);
743 }
745 }
746
747private:
748 boost::shared_ptr<PreviousData> commonData;
750 std::map<int, BlockData> setOfBlock;
751};
752
753// 10. Assembly of tangent for F_v_tau
754struct OpAssembleLhsVTau : OpFaceEle // C_VTau
755{
756 OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
757 : OpFaceEle(mass_field, flux_field, OpFaceEle::OPROWCOL) {
758 sYmm = false;
759 }
760
761 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
762 EntityType col_type, EntData &row_data,
763 EntData &col_data) {
765 const int nb_row_dofs = row_data.getIndices().size();
766 const int nb_col_dofs = col_data.getIndices().size();
767
768 if (nb_row_dofs && nb_col_dofs) {
769 mat.resize(nb_row_dofs, nb_col_dofs, false);
770 mat.clear();
771 const int nb_integration_pts = getGaussPts().size2();
772 auto t_w = getFTensor0IntegrationWeight();
773 auto t_row_v_base = row_data.getFTensor0N();
774 const double vol = getMeasure();
775
776 for (int gg = 0; gg != nb_integration_pts; ++gg) {
777 const double a = vol * t_w;
778 for (int rr = 0; rr != nb_row_dofs; ++rr) {
779 auto t_col_tau_grad = col_data.getFTensor2DiffN<3, 2>(gg, 0);
780 for (int cc = 0; cc != nb_col_dofs; ++cc) {
781 double div_col_base = t_col_tau_grad(0, 0) + t_col_tau_grad(1, 1);
782 mat(rr, cc) += (t_row_v_base * div_col_base) * a;
783 ++t_col_tau_grad;
784 }
785 ++t_row_v_base;
786 }
787 ++t_w;
788 }
789 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
790 ADD_VALUES);
791 }
793 }
794
795private:
797};
798
799// 11. Assembly of tangent for F_v_v
800struct OpAssembleLhsVV : OpFaceEle // D
801{
802 typedef boost::function<double(const double, const double)> FUval;
803 OpAssembleLhsVV(std::string mass_field,
804 boost::shared_ptr<PreviousData> &datau,
805 boost::shared_ptr<PreviousData> &datav,
806 FUval Drhs_u)
807 : OpFaceEle(mass_field, mass_field, OpFaceEle::OPROWCOL)
808 , DRhs_u(Drhs_u)
809 , dataU(datau)
810 , dataV(datav){
811 sYmm = true;
812 }
813
814 boost::shared_ptr<PreviousData> &dataU;
815 boost::shared_ptr<PreviousData> &dataV;
816
818
819 MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
820 EntityType col_type, EntData &row_data,
821 EntData &col_data) {
823
824 const int nb_row_dofs = row_data.getIndices().size();
825 const int nb_col_dofs = col_data.getIndices().size();
826 if (nb_row_dofs && nb_col_dofs) {
827
828 mat.resize(nb_row_dofs, nb_col_dofs, false);
829 mat.clear();
830 const int nb_integration_pts = getGaussPts().size2();
831
832 auto t_row_v_base = row_data.getFTensor0N();
833 auto t_u_value = getFTensor0FromVec(dataU->mass_values);
834 auto t_v_value = getFTensor0FromVec(dataV->mass_values);
835
836 auto t_w = getFTensor0IntegrationWeight();
837 const double ts_a = getFEMethod()->ts_a;
838 const double vol = getMeasure();
839
840 for (int gg = 0; gg != nb_integration_pts; ++gg) {
841 const double a = vol * t_w;
842 double dfu = DRhs_u(t_u_value, t_v_value);
843 for (int rr = 0; rr != nb_row_dofs; ++rr) {
844 auto t_col_v_base = col_data.getFTensor0N(gg, 0);
845 for (int cc = 0; cc != nb_col_dofs; ++cc) {
846 mat(rr, cc) += ((ts_a - 0*dfu) * t_row_v_base * t_col_v_base) * a;
847
848 ++t_col_v_base;
849 }
850 ++t_row_v_base;
851
852 }
853 ++t_w;
854 ++t_u_value;
855 ++t_v_value;
856 }
857 CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
858 ADD_VALUES);
859 if (row_side != col_side || row_type != col_type) {
860 transMat.resize(nb_col_dofs, nb_row_dofs, false);
861 noalias(transMat) = trans(mat);
862 CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
863 &transMat(0, 0), ADD_VALUES);
864 }
865 }
867 }
868
869private:
871};
872
873// struct OpError : public OpFaceEle {
874// typedef boost::function<double(const double, const double, const double)>
875// FVal;
876// typedef boost::function<FTensor::Tensor1<double, 3>(
877// const double, const double, const double)>
878// FGrad;
879// double &eRror;
880// OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad,
881// boost::shared_ptr<PreviousData> &prev_data,
882// std::map<int, BlockData> &block_map, double &err)
883// : OpFaceEle("ERROR", OpFaceEle::OPROW), exactVal(exact_value),
884// exactLap(exact_lap), exactGrad(exact_grad), prevData(prev_data),
885// setOfBlock(block_map), eRror(err) {}
886// MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
887// MoFEMFunctionBegin;
888// const int nb_dofs = data.getFieldData().size();
889// // cout << "nb_error_dofs : " << nb_dofs << endl;
890// if (nb_dofs) {
891// auto find_block_data = [&]() {
892// EntityHandle fe_ent = getFEEntityHandle();
893// BlockData *block_raw_ptr = nullptr;
894// for (auto &m : setOfBlock) {
895// if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end())
896// {
897// block_raw_ptr = &m.second;
898// break;
899// }
900// }
901// return block_raw_ptr;
902// };
903
904// auto block_data_ptr = find_block_data();
905// if (!block_data_ptr)
906// SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not
907// found");
908
909// auto &block_data = *block_data_ptr;
910
911// auto t_flux_value = getFTensor1FromMat<3>(prevData->flux_values);
912// // auto t_mass_dot = getFTensor0FromVec(prevData->mass_dots);
913// auto t_mass_value = getFTensor0FromVec(prevData->mass_values);
914// // auto t_flux_div = getFTensor0FromVec(prevData->flux_divs);
915// data.getFieldData().clear();
916// const double vol = getMeasure();
917// const int nb_integration_pts = getGaussPts().size2();
918// auto t_w = getFTensor0IntegrationWeight();
919// double dt;
920// CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
921// double ct = getFEMethod()->ts_t - dt;
922// auto t_coords = getFTensor1CoordsAtGaussPts();
923
924// FTensor::Tensor1<double, 3> t_exact_flux, t_flux_error;
925
926// for (int gg = 0; gg != nb_integration_pts; ++gg) {
927// const double a = vol * t_w;
928// double mass_exact = exactVal(t_coords(NX), t_coords(NY), ct);
929// // double flux_lap = - block_data.B0 * exactLap(t_coords(NX),
930// // t_coords(NY), ct);
931// t_exact_flux(i) =
932// -block_data.B0 * exactGrad(t_coords(NX), t_coords(NY), ct)(i);
933// t_flux_error(0) = t_flux_value(0) - t_exact_flux(0);
934// t_flux_error(1) = t_flux_value(1) - t_exact_flux(1);
935// t_flux_error(2) = t_flux_value(2) - t_exact_flux(2);
936// double local_error = pow(mass_exact - t_mass_value, 2) +
937// t_flux_error(i) * t_flux_error(i);
938// // cout << "flux_div : " << t_flux_div << " flux_exact : " <<
939// // flux_exact << endl;
940// data.getFieldData()[0] += a * local_error;
941// eRror += a * local_error;
942
943// ++t_w;
944// ++t_mass_value;
945// // ++t_flux_div;
946// ++t_flux_value;
947// // ++t_mass_dot;
948// ++t_coords;
949// }
950
951// data.getFieldDofs()[0]->getFieldData() = data.getFieldData()[0];
952// }
953// MoFEMFunctionReturn(0);
954// }
955
956// private:
957// FVal exactVal;
958// FVal exactLap;
959// FGrad exactGrad;
960// boost::shared_ptr<PreviousData> prevData;
961// std::map<int, BlockData> setOfBlock;
962
963// FTensor::Number<0> NX;
964// FTensor::Number<1> NY;
965// };
966
967struct Monitor : public FEMethod {
968 Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj<DM> &dm,
969 boost::shared_ptr<PostProc> &post_proc)
970 : cOmm(comm), rAnk(rank), dM(dm), postProc(post_proc){};
975 if (ts_step % save_every_nth_step == 0) {
977 CHKERR postProc->writeFile(
978 "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
979 }
980
982 }
983
984private:
985 SmartPetscObj<DM> dM;
986
987 boost::shared_ptr<PostProc> postProc;
988 MPI_Comm cOmm;
989 const int rAnk;
990};
991
992}; // namespace ElectroPhysiology
993
994#endif //__ELECPHYSOPERATORS_HPP__
static Index< 'K', 3 > K
ForcesAndSourcesCore::UserDataOperator UserDataOperator
constexpr double a
EntitiesFieldData::EntData EntData
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
size_t cholesky_decompose(const MATRIX &A, TRIA &L)
decompose the symmetric positive definit matrix A into product L L^T.
Definition: cholesky.hpp:52
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#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 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
double dt
Definition: heat_method.cpp:40
const double T
DataForcesAndSourcesCore::EntData EntData
const double essen_value
MoFEM::EdgeElementForcesAndSourcesCoreSwitch< 0 > BoundaryEle
FTensor::Index< 'i', 3 > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
const double r
rate factor
double duration
impulse duration (ns)
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProc > &post_proc)
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
boost::shared_ptr< PostProc > postProc
boost::shared_ptr< PreviousData > commonData
OpAssembleLhsTauTau(std::string flux_field, boost::shared_ptr< PreviousData > &commonData, std::map< int, BlockData > &block_map)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsTauV(std::string flux_field, std::string mass_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
boost::shared_ptr< PreviousData > commonData
OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsVV(std::string mass_field, boost::shared_ptr< PreviousData > &datau, boost::shared_ptr< PreviousData > &datav, FUval Drhs_u)
boost::function< double(const double, const double)> FUval
boost::shared_ptr< PreviousData > & dataV
boost::shared_ptr< PreviousData > & dataU
boost::shared_ptr< PreviousData > commonDatav
boost::function< double(const double, const double)> FUVal
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonDatau
OpAssembleSlowRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &common_datau, boost::shared_ptr< PreviousData > &common_datav, FUVal rhs_u)
OpAssembleStiffRhsTau(std::string flux_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonData
boost::function< double(const double, const double)> FUval
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonDatau
boost::shared_ptr< PreviousData > commonDatav
OpAssembleStiffRhsV(std::string flux_field, boost::shared_ptr< PreviousData > &datau, boost::shared_ptr< PreviousData > &datav, FUval rhs_u, std::map< int, BlockData > &block_map, Range &stim_region)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpEssentialBC(const std::string &flux_field, Range &essential_bd_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpInitialMass(const std::string &mass_field, Range &inner_surface, double &init_val)
OpSolveRecovery(const std::string &mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Method runge_kutta4)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< double(const double, const double, const double)> Method
boost::shared_ptr< PreviousData > dataV
boost::shared_ptr< PreviousData > dataU
bool sYmm
If true assume that matrix is symmetric structure.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
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
PetscReal ts_t
time
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Monitor solution.
Postprocess on face.