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