v0.10.0
ElecPhysOperators.hpp
Go to the documentation of this file.
1 #ifndef __ELECPHYSOPERATORS_HPP__
2 #define __ELECPHYSOPERATORS_HPP__
3 
4 #include <stdlib.h>
5 #include <BasicFiniteElements.hpp>
6 
7 namespace ElectroPhysiology {
8 
10 
12 
15 
17 
18 
19 
20 const int save_every_nth_step = 1;
21 
22 const double essen_value = 0;
23 
25 
26 
27 // problem parameters
28 const double alpha = 0.01;
29 const double gma = 0.002;
30 const double b = 0.15;
31 const double c = 8.00;
32 const double mu1 = 0.20;
33 const double mu2 = 0.30;
34 
35 const double D_tilde = 0.01;
36 
37 
38 
39 struct PreviousData {
42 
45 
47 
50 
52  jac.resize(3, 3, false);
53  inv_jac.resize(3, 3, false);
54  }
55 };
56 
57 
58 struct OpEssentialBC : public OpFaceEle {
59  OpEssentialBC(std::string flux_field, Range &essential_bd_ents)
61  }
62  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
64  int nb_dofs = data.getIndices().size();
65 
66  if (nb_dofs) {
68  bool is_essential =
69  (essential_bd_ents.find(fe_ent) != essential_bd_ents.end());
70  if (is_essential) {
71  int nb_gauss_pts = getGaussPts().size2();
72  int size2 = data.getN().size2();
73  if (3 * nb_dofs != static_cast<int>(data.getN().size2()))
74  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
75  "wrong number of dofs");
76  nN.resize(nb_dofs, nb_dofs, false);
77  nF.resize(nb_dofs, false);
78  nN.clear();
79  nF.clear();
80 
81  auto t_row_tau = data.getFTensor1N<3>();
82 
83  double *normal_ptr;
84  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
85  // HO geometry
86  normal_ptr = &getNormalsAtGaussPts(0)[0];
87  } else {
88  // Linear geometry, i.e. constant normal on face
89  normal_ptr = &getNormal()[0];
90  }
91  // set tensor from pointer
92  FTensor::Tensor1<const double *, 3> t_normal(normal_ptr, &normal_ptr[1],
93  &normal_ptr[2], 3);
94 
95  auto t_w = getFTensor0IntegrationWeight();
96  const double vol = getMeasure();
97  double nrm2 = 0;
98  for (int gg = 0; gg < nb_gauss_pts; gg++) {
99  if (gg == 0) {
100  nrm2 = sqrt(t_normal(i) * t_normal(i));
101  }
102  const double a = t_w * vol;
103  for (int rr = 0; rr != nb_dofs; rr++) {
105  &data.getVectorN<3>(gg)(0, HVEC0),
106  &data.getVectorN<3>(gg)(0, HVEC1),
107  &data.getVectorN<3>(gg)(0, HVEC2), 3);
108  nF[rr] += a * essen_value * t_row_tau(i) * t_normal(i) / nrm2;
109  for (int cc = 0; cc != nb_dofs; cc++) {
110  nN(rr, cc) += a * (t_row_tau(i) * t_normal(i)) *
111  (t_col_tau(i) * t_normal(i)) / (nrm2 * nrm2);
112  ++t_col_tau;
113  }
114  ++t_row_tau;
115  }
116  // If HO geometry increment t_normal to next integration point
117  if (getNormalsAtGaussPts().size1() == (unsigned int)nb_gauss_pts) {
118  ++t_normal;
119  nrm2 = sqrt(t_normal(i) * t_normal(i));
120  }
121  ++t_w;
122  }
123 
125  cholesky_solve(nN, nF, ublas::lower());
126 
127  for (auto &dof : data.getFieldDofs()) {
128  dof->getFieldData() = nF[dof->getEntDofIdx()];
129  }
130  }
131  }
133  }
134 
135 private:
139 };
140 
141 // Assembly of system mass matrix
142 // //***********************************************
143 
144 // Mass matrix corresponding to the flux equation.
145 // 01. Note that it is an identity matrix
146 
147 struct OpInitialMass : public OpVolEle {
148  OpInitialMass(const std::string &mass_field, Range &inner_surface,
149  double init_val)
150  : OpVolEle(mass_field, OpVolEle::OPROW), innerSurface(inner_surface),
151  initVal(init_val) {}
154  Range &innerSurface;
155  double initVal;
156  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
158  int nb_dofs = data.getFieldData().size();
159  if (nb_dofs) {
160  EntityHandle fe_ent = getFEEntityHandle();
161  bool is_inner_side = (innerSurface.find(fe_ent) != innerSurface.end());
162  if (is_inner_side) {
163  int nb_gauss_pts = getGaussPts().size2();
164  if (nb_dofs != static_cast<int>(data.getN().size2()))
165  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
166  "wrong number of dofs");
167  nN.resize(nb_dofs, nb_dofs, false);
168  nF.resize(nb_dofs, false);
169  nN.clear();
170  nF.clear();
171 
172  auto t_row_mass = data.getFTensor0N();
173  auto t_w = getFTensor0IntegrationWeight();
174  const double vol = getMeasure();
175 
176  for (int gg = 0; gg < nb_gauss_pts; gg++) {
177  const double a = t_w * vol;
178  // double r = ((double)rand() / (RAND_MAX));
179  for (int rr = 0; rr != nb_dofs; rr++) {
180  auto t_col_mass = data.getFTensor0N(gg, 0);
181  nF[rr] += a * initVal * t_row_mass;
182  for (int cc = 0; cc != nb_dofs; cc++) {
183  nN(rr, cc) += a * t_row_mass * t_col_mass;
184  ++t_col_mass;
185  }
186  ++t_row_mass;
187  }
188  ++t_w;
189  }
190 
192  cholesky_solve(nN, nF, ublas::lower());
193 
194  for (auto &dof : data.getFieldDofs()) {
195  dof->getFieldData() = nF[dof->getEntDofIdx()];
196 
197  // this is only to check
198  // data.getFieldData()[dof->getEntDofIdx()] = nF[dof->getEntDofIdx()];
199  }
200  }
201  }
203  }
204 };
205 
206 
207 struct OpSolveRecovery : public OpVolEle {
208  typedef boost::function<double(const double, const double, const double)>
210  OpSolveRecovery(const std::string &mass_field,
211  boost::shared_ptr<PreviousData> &data_u,
212  boost::shared_ptr<PreviousData> &data_v,
213  Method runge_kutta4)
214  : OpVolEle(mass_field, OpVolEle::OPROW)
215  , dataU(data_u)
216  , dataV(data_v)
217  , rungeKutta4(runge_kutta4)
218  {}
219  boost::shared_ptr<PreviousData> dataU;
220  boost::shared_ptr<PreviousData> dataV;
222 
225  double initVal;
226  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
228  int nb_dofs = data.getFieldData().size();
229  if (nb_dofs) {
230  int nb_gauss_pts = getGaussPts().size2();
231  if (nb_dofs != static_cast<int>(data.getN().size2()))
232  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
233  "wrong number of dofs");
234  nN.resize(nb_dofs, nb_dofs, false);
235  nF.resize(nb_dofs, false);
236  nN.clear();
237  nF.clear();
238 
239  auto t_val_u = getFTensor0FromVec(dataU->mass_values);
240  auto t_val_v = getFTensor0FromVec(dataV->mass_values);
241 
242  double dt;
243  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
244 
245  auto t_row_mass = data.getFTensor0N();
246  auto t_w = getFTensor0IntegrationWeight();
247  const double vol = getMeasure();
248 
249  for (int gg = 0; gg < nb_gauss_pts; gg++) {
250  const double a = t_w * vol;
251  const double vn = rungeKutta4(t_val_u, t_val_v, dt);
252  for (int rr = 0; rr != nb_dofs; rr++) {
253  auto t_col_mass = data.getFTensor0N(gg, 0);
254  nF[rr] += a * vn * t_row_mass;
255  for (int cc = 0; cc != nb_dofs; cc++) {
256  nN(rr, cc) += a * t_row_mass * t_col_mass;
257  ++t_col_mass;
258  }
259  ++t_row_mass;
260  }
261  ++t_w;
262  ++t_val_u;
263  ++t_val_v;
264  }
265 
267  cholesky_solve(nN, nF, ublas::lower());
268 
269  for (auto &dof : data.getFieldDofs()) {
270  dof->getFieldData() = nF[dof->getEntDofIdx()];
271 
272  // this is only to check
273  // data.getFieldData()[dof->getEntDofIdx()] = nF[dof->getEntDofIdx()];
274  }
275 
276  }
278  }
279 };
280 
281 // Assembly of RHS for explicit (slow)
282 // part//**************************************
283 
284 // 2. RHS for explicit part of the mass balance equation
286 {
287  typedef boost::function<double(const double, const double)>
289  OpAssembleSlowRhsV(std::string mass_field,
290  boost::shared_ptr<PreviousData> &data_u,
291  boost::shared_ptr<PreviousData> &data_v,
292  Feval_u rhs_u)
293  : OpVolEle(mass_field, OpVolEle::OPROW)
294  , dataU(data_u)
295  , dataV(data_v)
296  , rhsU(rhs_u) {}
297 
302 
303  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
305 
306  const int nb_dofs = data.getIndices().size();
307  if (nb_dofs) {
308  if (nb_dofs != static_cast<int>(data.getN().size2()))
309  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
310  "wrong number of dofs");
311  vecF.resize(nb_dofs, false);
312  mat.resize(nb_dofs, nb_dofs, false);
313  vecF.clear();
314  mat.clear();
315  const int nb_integration_pts = getGaussPts().size2();
316  auto t_val_u = getFTensor0FromVec(dataU->mass_values);
317  auto t_val_v = getFTensor0FromVec(dataV->mass_values);
318 
319 
320  auto t_row_v_base = data.getFTensor0N();
321  auto t_w = getFTensor0IntegrationWeight();
322  const double vol = getMeasure();
323 
324 
325  for (int gg = 0; gg != nb_integration_pts; ++gg) {
326 
327  double rhs = rhsU(t_val_u, t_val_v);
328  const double a = vol * t_w;
329 
330  for (int rr = 0; rr != nb_dofs; ++rr) {
331  auto t_col_v_base = data.getFTensor0N(gg, 0);
332  vecF[rr] += a * rhs * t_row_v_base;
333  for (int cc = 0; cc != nb_dofs; ++cc) {
334  mat(rr, cc) += a * t_row_v_base * t_col_v_base;
335  ++t_col_v_base;
336  }
337  ++t_row_v_base;
338  }
339  ++t_val_u;
340  ++t_val_v;
341  ++t_w;
342 
343  }
345  cholesky_solve(mat, vecF, ublas::lower());
346 
347  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
348  PETSC_TRUE);
349  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
350  ADD_VALUES);
351  }
353  }
354 
355 private:
356  boost::shared_ptr<PreviousData> dataU;
357  boost::shared_ptr<PreviousData> dataV;
360 
361 };
362 
363 // // 5. RHS contribution of the natural boundary condition
364 // struct OpAssembleNaturalBCRhsTau : OpFaceEle // R_tau_2
365 // {
366 // OpAssembleNaturalBCRhsTau(std::string flux_field, Range &natural_bd_ents)
367 // : OpFaceEle(flux_field, OpFaceEle::OPROW),
368 // natural_bd_ents(natural_bd_ents) {}
369 
370 // MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
371 // MoFEMFunctionBegin;
372 // const int nb_dofs = data.getIndices().size();
373 
374 // if (nb_dofs) {
375 // EntityHandle row_side_ent = data.getFieldDofs()[0]->getEnt();
376 
377 // bool is_natural =
378 // (natural_bd_ents.find(row_side_ent) != natural_bd_ents.end());
379 // if (is_natural) {
380 // // cerr << "In NaturalBCRhsTau..." << endl;
381 // vecF.resize(nb_dofs, false);
382 // vecF.clear();
383 // const int nb_integration_pts = getGaussPts().size2();
384 // auto t_tau_base = data.getFTensor1N<3>();
385 
386 // auto dir = getDirection();
387 // FTensor::Tensor1<double, 3> t_normal(-dir[1], dir[0], dir[2]);
388 
389 // auto t_w = getFTensor0IntegrationWeight();
390 
391 // for (int gg = 0; gg != nb_integration_pts; ++gg) {
392 // const double a = t_w;
393 // for (int rr = 0; rr != nb_dofs; ++rr) {
394 // vecF[rr] += (t_tau_base(i) * t_normal(i) * natu_value) * a;
395 // ++t_tau_base;
396 // }
397 // ++t_w;
398 // }
399 // CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
400 // PETSC_TRUE);
401 // CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
402 // ADD_VALUES);
403 // }
404 // }
405 // MoFEMFunctionReturn(0);
406 // }
407 
408 // private:
409 // VectorDouble vecF;
410 // Range natural_bd_ents;
411 // };
412 
413 // Assembly of RHS for the implicit (stiff) part excluding the essential
414 // boundary //**********************************
415 // 3. Assembly of F_tau excluding the essential boundary condition
416 template <int dim>
417 struct OpAssembleStiffRhsTau : OpVolEle // F_tau_1
418 {
419  OpAssembleStiffRhsTau(std::string flux_field,
420  boost::shared_ptr<PreviousData> &data)
421  : OpVolEle(flux_field, OpVolEle::OPROW), commonData(data) {}
422 
423  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
425 
426  const int nb_dofs = data.getIndices().size();
427  if (nb_dofs) {
428 
429 
430  vecF.resize(nb_dofs, false);
431  vecF.clear();
432 
433  const int nb_integration_pts = getGaussPts().size2();
434  auto t_flux_value = getFTensor1FromMat<3>(commonData->flux_values);
435  auto t_mass_value = getFTensor0FromVec(commonData->mass_values);
436  auto t_tau_base = data.getFTensor1N<3>();
437 
438  auto t_tau_grad = data.getFTensor2DiffN<3, 3>();
439 
440  auto t_w = getFTensor0IntegrationWeight();
441  const double vol = getMeasure();
442 
443  for (int gg = 0; gg < nb_integration_pts; ++gg) {
444 
445  const double K_inv = 1. / D_tilde;
446  const double a = vol * t_w;
447  for (int rr = 0; rr < nb_dofs; ++rr) {
448  double div_base =
449  t_tau_grad(0, 0) + t_tau_grad(1, 1) + t_tau_grad(2, 2);
450  vecF[rr] += (K_inv * t_tau_base(i) * t_flux_value(i) -
451  div_base * t_mass_value) *
452  a;
453  ++t_tau_base;
454  ++t_tau_grad;
455  }
456  ++t_flux_value;
457  ++t_mass_value;
458  ++t_w;
459  }
460  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
461  PETSC_TRUE);
462  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
463  ADD_VALUES);
464  }
466  }
467 
468 private:
469  boost::shared_ptr<PreviousData> commonData;
471 };
472 // 4. Assembly of F_v
473 template <int dim>
475 {
476  OpAssembleStiffRhsV(std::string mass_field,
477  boost::shared_ptr<PreviousData> &data_u,
478  Range &stim_region)
479  : OpVolEle(mass_field, OpVolEle::OPROW)
480  , dataU(data_u)
481  , stimRegion(stim_region) {}
482 
483  Range &stimRegion;
484 
485  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
487  const int nb_dofs = data.getIndices().size();
488  // cerr << "In StiffRhsV ..." << endl;
489  if (nb_dofs) {
490 
491  vecF.resize(nb_dofs, false);
492  vecF.clear();
493  const int nb_integration_pts = getGaussPts().size2();
494  auto t_dot_u = getFTensor0FromVec(dataU->mass_dots);
495 
496 
497  auto t_div_u = getFTensor0FromVec(dataU->flux_divs);
498  auto t_row_v_base = data.getFTensor0N();
499  auto t_w = getFTensor0IntegrationWeight();
500  const double vol = getMeasure();
501  const double c_time = getFEMethod()->ts_t;
502 
503  double dt;
504  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
505 
506  double stim = 0.0;
507 
508  double T = 627.0;
509  double duration = 5.0;
510 
511  if (T - dt < c_time && c_time <= T + duration) {
512  EntityHandle stim_ent = getFEEntityHandle();
513  if (stimRegion.find(stim_ent) != stimRegion.end()) {
514  stim = 40.0;
515  } else {
516  stim = 0.0;
517  }
518  }
519 
520  for (int gg = 0; gg < nb_integration_pts; ++gg) {
521  const double a = vol * t_w;
522  for (int rr = 0; rr < nb_dofs; ++rr) {
523  vecF[rr] += t_row_v_base * (t_dot_u + t_div_u - stim) * a;
524  ++t_row_v_base;
525  }
526  ++t_dot_u;
527  ++t_div_u;
528  ++t_w;
529  }
530  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
531  PETSC_TRUE);
532  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
533  ADD_VALUES);
534  }
536  }
537 
538 private:
539  boost::shared_ptr<PreviousData> dataU;
540  boost::shared_ptr<PreviousData> dataV;
542 };
543 
544 // Tangent operator
545 // //**********************************************
546 // 7. Tangent assembly for F_tautau excluding the essential boundary condition
547 template <int dim>
548 struct OpAssembleLhsTauTau : OpVolEle // A_TauTau_1
549 {
550  OpAssembleLhsTauTau(std::string flux_field)
551  : OpVolEle(flux_field, flux_field, OpVolEle::OPROWCOL) {
552  sYmm = true;
553  }
554 
555  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
556  EntityType col_type, EntData &row_data,
557  EntData &col_data) {
559  const int nb_row_dofs = row_data.getIndices().size();
560  const int nb_col_dofs = col_data.getIndices().size();
561 
562  if (nb_row_dofs && nb_col_dofs) {
563 
564  mat.resize(nb_row_dofs, nb_col_dofs, false);
565  mat.clear();
566  const int nb_integration_pts = getGaussPts().size2();
567 
568  auto t_row_tau_base = row_data.getFTensor1N<3>();
569 
570  auto t_w = getFTensor0IntegrationWeight();
571  const double vol = getMeasure();
572 
573  for (int gg = 0; gg != nb_integration_pts; ++gg) {
574  const double a = vol * t_w;
575  const double K_inv = 1. / D_tilde;
576  for (int rr = 0; rr != nb_row_dofs; ++rr) {
577  auto t_col_tau_base = col_data.getFTensor1N<3>(gg, 0);
578  for (int cc = 0; cc != nb_col_dofs; ++cc) {
579  mat(rr, cc) += (K_inv * t_row_tau_base(i) * t_col_tau_base(i)) * a;
580  ++t_col_tau_base;
581  }
582  ++t_row_tau_base;
583  }
584  ++t_w;
585  }
586  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
587  ADD_VALUES);
588  if (row_side != col_side || row_type != col_type) {
589  transMat.resize(nb_col_dofs, nb_row_dofs, false);
590  noalias(transMat) = trans(mat);
591  CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
592  &transMat(0, 0), ADD_VALUES);
593  }
594  }
596  }
597 
598 private:
600 };
601 
602 // 9. Assembly of tangent for F_tau_v excluding the essential bc
603 template <int dim>
604 struct OpAssembleLhsTauV : OpVolEle // E_TauV
605 {
606  OpAssembleLhsTauV(std::string flux_field, std::string mass_field)
607  : OpVolEle(flux_field, mass_field, OpVolEle::OPROWCOL) {
608  sYmm = false;
609  }
610 
611  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
612  EntityType col_type, EntData &row_data,
613  EntData &col_data) {
615  const int nb_row_dofs = row_data.getIndices().size();
616  const int nb_col_dofs = col_data.getIndices().size();
617 
618  if (nb_row_dofs && nb_col_dofs) {
619 
620  mat.resize(nb_row_dofs, nb_col_dofs, false);
621  mat.clear();
622  const int nb_integration_pts = getGaussPts().size2();
623  auto t_w = getFTensor0IntegrationWeight();
624  auto t_row_tau_base = row_data.getFTensor1N<3>();
625 
626  auto t_row_tau_grad = row_data.getFTensor2DiffN<3, 3>();
627  const double vol = getMeasure();
628 
629  for (int gg = 0; gg != nb_integration_pts; ++gg) {
630  const double a = vol * t_w;
631 
632  for (int rr = 0; rr != nb_row_dofs; ++rr) {
633  auto t_col_v_base = col_data.getFTensor0N(gg, 0);
634  for (int cc = 0; cc != nb_col_dofs; ++cc) {
635  double div_row_base = t_row_tau_grad(0, 0) + t_row_tau_grad(1, 1) +
636  t_row_tau_grad(2, 2);
637  mat(rr, cc) += -(div_row_base * t_col_v_base) * a;
638  ++t_col_v_base;
639  }
640  ++t_row_tau_base;
641  ++t_row_tau_grad;
642  }
643  ++t_w;
644  }
645  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
646  ADD_VALUES);
647  }
649  }
650 
651 private:
653 };
654 
655 // 10. Assembly of tangent for F_v_tau
656 struct OpAssembleLhsVTau : OpVolEle // C_VTau
657 {
658  OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
659  : OpVolEle(mass_field, flux_field, OpVolEle::OPROWCOL) {
660  sYmm = false;
661  }
662 
663  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
664  EntityType col_type, EntData &row_data,
665  EntData &col_data) {
667  const int nb_row_dofs = row_data.getIndices().size();
668  const int nb_col_dofs = col_data.getIndices().size();
669 
670  if (nb_row_dofs && nb_col_dofs) {
671  mat.resize(nb_row_dofs, nb_col_dofs, false);
672  mat.clear();
673  const int nb_integration_pts = getGaussPts().size2();
674  auto t_w = getFTensor0IntegrationWeight();
675  auto t_row_v_base = row_data.getFTensor0N();
676  const double vol = getMeasure();
677 
678  for (int gg = 0; gg != nb_integration_pts; ++gg) {
679  const double a = vol * t_w;
680  for (int rr = 0; rr != nb_row_dofs; ++rr) {
681  auto t_col_tau_grad = col_data.getFTensor2DiffN<3, 3>(gg, 0);
682  for (int cc = 0; cc != nb_col_dofs; ++cc) {
683  double div_col_base = t_col_tau_grad(0, 0) + t_col_tau_grad(1, 1) + t_col_tau_grad(2, 2);
684  mat(rr, cc) += (t_row_v_base * div_col_base) * a;
685  ++t_col_tau_grad;
686  }
687  ++t_row_v_base;
688  }
689  ++t_w;
690  }
691  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
692  ADD_VALUES);
693  }
695  }
696 
697 private:
699 };
700 
701 // 11. Assembly of tangent for F_v_v
703 {
704  OpAssembleLhsVV(std::string mass_field)
705  : OpVolEle(mass_field, mass_field, OpVolEle::OPROWCOL) {
706  sYmm = true;
707  }
708 
709  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
710  EntityType col_type, EntData &row_data,
711  EntData &col_data) {
713 
714  const int nb_row_dofs = row_data.getIndices().size();
715  const int nb_col_dofs = col_data.getIndices().size();
716  if (nb_row_dofs && nb_col_dofs) {
717 
718  mat.resize(nb_row_dofs, nb_col_dofs, false);
719  mat.clear();
720  const int nb_integration_pts = getGaussPts().size2();
721 
722  auto t_row_v_base = row_data.getFTensor0N();
723 
724  auto t_w = getFTensor0IntegrationWeight();
725  const double ts_a = getFEMethod()->ts_a;
726  const double vol = getMeasure();
727 
728  for (int gg = 0; gg != nb_integration_pts; ++gg) {
729  const double a = vol * t_w;
730 
731  for (int rr = 0; rr != nb_row_dofs; ++rr) {
732  auto t_col_v_base = col_data.getFTensor0N(gg, 0);
733  for (int cc = 0; cc != nb_col_dofs; ++cc) {
734  mat(rr, cc) += (ts_a * t_row_v_base * t_col_v_base) * a;
735 
736  ++t_col_v_base;
737  }
738  ++t_row_v_base;
739  }
740  ++t_w;
741  }
742  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
743  ADD_VALUES);
744  if (row_side != col_side || row_type != col_type) {
745  transMat.resize(nb_col_dofs, nb_row_dofs, false);
746  noalias(transMat) = trans(mat);
747  CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
748  &transMat(0, 0), ADD_VALUES);
749  }
750  }
752  }
753 
754 private:
756 };
757 
758 // struct OpError : public OpFaceEle {
759 // typedef boost::function<double(const double, const double, const double)>
760 // FVal;
761 // typedef boost::function<FTensor::Tensor1<double, 3>(
762 // const double, const double, const double)>
763 // FGrad;
764 // double &eRror;
765 // OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad,
766 // boost::shared_ptr<PreviousData> &prev_data,
767 // std::map<int, BlockData> &block_map, double &err)
768 // : OpFaceEle("ERROR", OpFaceEle::OPROW), exactVal(exact_value),
769 // exactLap(exact_lap), exactGrad(exact_grad), prevData(prev_data),
770 // setOfBlock(block_map), eRror(err) {}
771 // MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
772 // MoFEMFunctionBegin;
773 // const int nb_dofs = data.getFieldData().size();
774 // // cout << "nb_error_dofs : " << nb_dofs << endl;
775 // if (nb_dofs) {
776 // auto find_block_data = [&]() {
777 // EntityHandle fe_ent = getFEEntityHandle();
778 // BlockData *block_raw_ptr = nullptr;
779 // for (auto &m : setOfBlock) {
780 // if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
781 // block_raw_ptr = &m.second;
782 // break;
783 // }
784 // }
785 // return block_raw_ptr;
786 // };
787 
788 // auto block_data_ptr = find_block_data();
789 // if (!block_data_ptr)
790 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
791 
792 // auto &block_data = *block_data_ptr;
793 
794 // auto t_flux_value = getFTensor1FromMat<3>(prevData->flux_values);
795 // // auto t_mass_dot = getFTensor0FromVec(prevData->mass_dots);
796 // auto t_mass_value = getFTensor0FromVec(prevData->mass_values);
797 // // auto t_flux_div = getFTensor0FromVec(prevData->flux_divs);
798 // data.getFieldData().clear();
799 // const double vol = getMeasure();
800 // const int nb_integration_pts = getGaussPts().size2();
801 // auto t_w = getFTensor0IntegrationWeight();
802 // double dt;
803 // CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
804 // double ct = getFEMethod()->ts_t - dt;
805 // auto t_coords = getFTensor1CoordsAtGaussPts();
806 
807 // FTensor::Tensor1<double, 3> t_exact_flux, t_flux_error;
808 
809 // for (int gg = 0; gg != nb_integration_pts; ++gg) {
810 // const double a = vol * t_w;
811 // double mass_exact = exactVal(t_coords(NX), t_coords(NY), ct);
812 // // double flux_lap = - block_data.B0 * exactLap(t_coords(NX),
813 // // t_coords(NY), ct);
814 // t_exact_flux(i) =
815 // -block_data.B0 * exactGrad(t_coords(NX), t_coords(NY), ct)(i);
816 // t_flux_error(0) = t_flux_value(0) - t_exact_flux(0);
817 // t_flux_error(1) = t_flux_value(1) - t_exact_flux(1);
818 // t_flux_error(2) = t_flux_value(2) - t_exact_flux(2);
819 // double local_error = pow(mass_exact - t_mass_value, 2) +
820 // t_flux_error(i) * t_flux_error(i);
821 // // cout << "flux_div : " << t_flux_div << " flux_exact : " <<
822 // // flux_exact << endl;
823 // data.getFieldData()[0] += a * local_error;
824 // eRror += a * local_error;
825 
826 // ++t_w;
827 // ++t_mass_value;
828 // // ++t_flux_div;
829 // ++t_flux_value;
830 // // ++t_mass_dot;
831 // ++t_coords;
832 // }
833 
834 // data.getFieldDofs()[0]->getFieldData() = data.getFieldData()[0];
835 // }
836 // MoFEMFunctionReturn(0);
837 // }
838 
839 // private:
840 // FVal exactVal;
841 // FVal exactLap;
842 // FGrad exactGrad;
843 // boost::shared_ptr<PreviousData> prevData;
844 // std::map<int, BlockData> setOfBlock;
845 
846 // FTensor::Number<0> NX;
847 // FTensor::Number<1> NY;
848 // };
849 
850 
851 
852 struct Monitor : public FEMethod {
853  Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj<DM> &dm,
854  boost::shared_ptr<PostProcVolumeOnRefinedMesh> &post_proc)
855  : cOmm(comm), rAnk(rank), dM(dm), postProc(post_proc){};
856  MoFEMErrorCode preProcess() { return 0; }
857  MoFEMErrorCode operator()() { return 0; }
860  if (ts_step % save_every_nth_step == 0) {
862  CHKERR postProc->writeFile(
863  "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
864  }
865 
867  }
868 
869 private:
870  SmartPetscObj<DM> dM;
871 
872  boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProc;
873  MPI_Comm cOmm;
874  const int rAnk;
875 };
876 
877 }; // namespace ReactionDiffusion
878 
879 #endif //__ELECPHYSOPERATORS_HPP__
ElectroPhysiology::OpAssembleStiffRhsV::dataV
boost::shared_ptr< PreviousData > dataV
Definition: ElecPhysOperators.hpp:540
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::OpAssembleStiffRhsV::OpAssembleStiffRhsV
OpAssembleStiffRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &data_u, Range &stim_region)
Definition: ElecPhysOperators.hpp:476
ElectroPhysiology::OpAssembleSlowRhsV::NZ
FTensor::Number< 2 > NZ
Definition: ElecPhysOperators.hpp:301
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
MoFEM::VolumeElementForcesAndSourcesCoreSwitch::UserDataOperator
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:344
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::PreviousData
Definition: ElecPhysOperators.hpp:39
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEM::Types::MatrixDouble
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
MoFEM::VolumeElementForcesAndSourcesCore
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
Definition: VolumeElementForcesAndSourcesCore.hpp:354
ElectroPhysiology::Monitor::operator()
MoFEMErrorCode operator()()
Definition: ElecPhysOperators.hpp:857
ElectroPhysiology::OpAssembleLhsTauV::OpAssembleLhsTauV
OpAssembleLhsTauV(std::string flux_field, std::string mass_field)
Definition: ElecPhysOperators.hpp:606
MoFEM::VolumeElementForcesAndSourcesCoreSwitch
Volume finite element with switches.
Definition: VolumeElementForcesAndSourcesCore.hpp:339
ElectroPhysiology::OpSolveRecovery::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: ElecPhysOperators.hpp:226
ElectroPhysiology::OpAssembleStiffRhsTau::vecF
VectorDouble vecF
Definition: ElecPhysOperators.hpp:470
FTensor::Tensor1
Definition: Tensor1_value.hpp:9
ElectroPhysiology::OpInitialMass
Definition: ElecPhysOperators.hpp:147
ElectroPhysiology::OpAssembleLhsTauTau::OpAssembleLhsTauTau
OpAssembleLhsTauTau(std::string flux_field)
Definition: ElecPhysOperators.hpp:550
ElectroPhysiology::OpAssembleSlowRhsV::NY
FTensor::Number< 1 > NY
Definition: ElecPhysOperators.hpp:300
ElectroPhysiology::OpSolveRecovery::initVal
double initVal
Definition: ElecPhysOperators.hpp:225
ElectroPhysiology::OpEssentialBC::OpEssentialBC
OpEssentialBC(std::string flux_field, Range &essential_bd_ents)
Definition: ElecPhysOperators.hpp:59
ElectroPhysiology::OpEssentialBC
Definition: ElecPhysOperators.hpp:58
ElectroPhysiology::Monitor::postProc
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
Definition: ElecPhysOperators.hpp:872
HVEC1
@ HVEC1
Definition: definitions.h:255
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
ElectroPhysiology::OpAssembleStiffRhsV::dataU
boost::shared_ptr< PreviousData > dataU
Definition: ElecPhysOperators.hpp:539
ElectroPhysiology::OpAssembleSlowRhsV::mat
MatrixDouble mat
Definition: ElecPhysOperators.hpp:359
ElectroPhysiology::PreviousData::mass_values
VectorDouble mass_values
Definition: ElecPhysOperators.hpp:44
ElectroPhysiology::OpInitialMass::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: ElecPhysOperators.hpp:156
ElectroPhysiology::OpAssembleStiffRhsTau::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: ElecPhysOperators.hpp:423
ElectroPhysiology::Monitor
Definition: ElecPhysOperators.hpp:852
ElectroPhysiology::OpEssentialBC::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: ElecPhysOperators.hpp:62
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::OpSolveRecovery::rungeKutta4
Method rungeKutta4
Definition: ElecPhysOperators.hpp:221
ElectroPhysiology::PreviousData::PreviousData
PreviousData()
Definition: ElecPhysOperators.hpp:51
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::FaceEle
MoFEM::FaceElementForcesAndSourcesCore FaceEle
Definition: ElecPhysOperators.hpp:11
ElectroPhysiology::alpha
const double alpha
Definition: ElecPhysOperators.hpp:28
ElectroPhysiology::Monitor::postProcess
MoFEMErrorCode postProcess()
Definition: ElecPhysOperators.hpp:858
T
ElectroPhysiology::OpInitialMass::OpInitialMass
OpInitialMass(const std::string &mass_field, Range &inner_surface, double init_val)
Definition: ElecPhysOperators.hpp:148
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: ElecPhysOperators.hpp:210
ElectroPhysiology::OpAssembleLhsTauV
Definition: ElecPhysOperators.hpp:605
FTensor::Index< 'i', 3 >
ElectroPhysiology::OpSolveRecovery::Method
boost::function< double(const double, const double, const double)> Method
Definition: ElecPhysOperators.hpp:209
ElectroPhysiology::PreviousData::slow_values
VectorDouble slow_values
Definition: ElecPhysOperators.hpp:46
ElectroPhysiology::b
const double b
Definition: ElecPhysOperators.hpp:30
ElectroPhysiology::OpAssembleSlowRhsV::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: ElecPhysOperators.hpp:303
MoFEM::TSMethod::ts_a
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
Definition: LoopMethods.hpp:198
ElectroPhysiology::OpAssembleLhsVV
Definition: ElecPhysOperators.hpp:703
ElectroPhysiology::D_tilde
const double D_tilde
Definition: ElecPhysOperators.hpp:35
ElectroPhysiology::OpAssembleLhsTauTau::transMat
MatrixDouble transMat
Definition: ElecPhysOperators.hpp:599
ElectroPhysiology::OpInitialMass::initVal
double initVal
Definition: ElecPhysOperators.hpp:155
ElectroPhysiology::OpAssembleStiffRhsTau::OpAssembleStiffRhsTau
OpAssembleStiffRhsTau(std::string flux_field, boost::shared_ptr< PreviousData > &data)
Definition: ElecPhysOperators.hpp:419
ElectroPhysiology::OpAssembleSlowRhsV::NX
FTensor::Number< 0 > NX
Definition: ElecPhysOperators.hpp:299
ElectroPhysiology::OpEssentialBC::nF
VectorDouble nF
Definition: ElecPhysOperators.hpp:137
MoFEM::ForcesAndSourcesCore::UserDataOperator::OPROWCOL
@ OPROWCOL
Definition: ForcesAndSourcesCore.hpp:101
ElectroPhysiology::Monitor::Monitor
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcVolumeOnRefinedMesh > &post_proc)
Definition: ElecPhysOperators.hpp:853
ElectroPhysiology::OpAssembleSlowRhsV::dataV
boost::shared_ptr< PreviousData > dataV
Definition: ElecPhysOperators.hpp:357
ElectroPhysiology::EntData
DataForcesAndSourcesCore::EntData EntData
Definition: ElecPhysOperators.hpp:16
MoFEM::ForcesAndSourcesCore::UserDataOperator::getFEMethod
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
Definition: ForcesAndSourcesCore.hpp:902
ElectroPhysiology::OpAssembleSlowRhsV::dataU
boost::shared_ptr< PreviousData > dataU
Definition: ElecPhysOperators.hpp:356
ElectroPhysiology::save_every_nth_step
const int save_every_nth_step
Definition: ElecPhysOperators.hpp:20
ElectroPhysiology::OpAssembleLhsVV::mat
MatrixDouble mat
Definition: ElecPhysOperators.hpp:755
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
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::Monitor::dM
SmartPetscObj< DM > dM
Definition: ElecPhysOperators.hpp:870
ElectroPhysiology::OpAssembleSlowRhsV
Definition: ElecPhysOperators.hpp:286
ElectroPhysiology::OpSolveRecovery::dataV
boost::shared_ptr< PreviousData > dataV
Definition: ElecPhysOperators.hpp:220
ElectroPhysiology::OpAssembleSlowRhsV::OpAssembleSlowRhsV
OpAssembleSlowRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Feval_u rhs_u)
Definition: ElecPhysOperators.hpp:289
ElectroPhysiology::OpAssembleLhsVTau
Definition: ElecPhysOperators.hpp:657
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: ElecPhysOperators.hpp:663
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
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::OpAssembleStiffRhsV
Definition: ElecPhysOperators.hpp:475
dt
constexpr double dt
Definition: heat_method.cpp:32
ElectroPhysiology::c
const double c
Definition: ElecPhysOperators.hpp:31
MoFEM::FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormalsAtGaussPts
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
Definition: FaceElementForcesAndSourcesCore.hpp:508
MoFEM::FaceElementForcesAndSourcesCore
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
Definition: FaceElementForcesAndSourcesCore.hpp:337
ElectroPhysiology::OpAssembleSlowRhsV::Feval_u
boost::function< double(const double, const double)> Feval_u
Definition: ElecPhysOperators.hpp:288
ElectroPhysiology::OpAssembleLhsTauTau::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Definition: ElecPhysOperators.hpp:555
ElectroPhysiology::PreviousData::flux_values
MatrixDouble flux_values
Definition: ElecPhysOperators.hpp:40
ElectroPhysiology::OpAssembleStiffRhsV::vecF
VectorDouble vecF
Definition: ElecPhysOperators.hpp:541
ElectroPhysiology::OpAssembleStiffRhsV::stimRegion
Range & stimRegion
Definition: ElecPhysOperators.hpp:483
HVEC2
@ HVEC2
Definition: definitions.h:255
MoFEM::VolumeElementForcesAndSourcesCoreBase::UserDataOperator
default operator for TET element
Definition: VolumeElementForcesAndSourcesCore.hpp:45
ElectroPhysiology::OpSolveRecovery
Definition: ElecPhysOperators.hpp:207
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: ElecPhysOperators.hpp:709
ElectroPhysiology::OpAssembleLhsVV::OpAssembleLhsVV
OpAssembleLhsVV(std::string mass_field)
Definition: ElecPhysOperators.hpp:704
MoFEM::FaceElementForcesAndSourcesCoreBase::UserDataOperator::getNormal
VectorDouble & getNormal()
get triangle normal
Definition: FaceElementForcesAndSourcesCore.hpp:424
ElectroPhysiology::OpAssembleSlowRhsV::rhsU
Feval_u rhsU
Definition: ElecPhysOperators.hpp:298
ElectroPhysiology::OpAssembleLhsTauV::doWork
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Definition: ElecPhysOperators.hpp:611
MoFEM::DataOperator::sYmm
bool sYmm
If true assume that matrix is symmetric structure.
Definition: DataOperators.hpp:72
convert.int
int
Definition: convert.py:66
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: ElecPhysOperators.hpp:658
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
ElectroPhysiology::OpAssembleLhsTauTau::mat
MatrixDouble mat
Definition: ElecPhysOperators.hpp:599
HVEC0
@ HVEC0
Definition: definitions.h:255
ElectroPhysiology::OpAssembleStiffRhsV::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Definition: ElecPhysOperators.hpp:485
ElectroPhysiology::Monitor::preProcess
MoFEMErrorCode preProcess()
Definition: ElecPhysOperators.hpp:856
MoFEM::ForcesAndSourcesCore::UserDataOperator::getGaussPts
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Definition: ForcesAndSourcesCore.hpp:1016
ElectroPhysiology::OpAssembleLhsVV::transMat
MatrixDouble transMat
Definition: ElecPhysOperators.hpp:755