v0.9.1
ElecPhysOperators2D.hpp
Go to the documentation of this file.
1 #ifndef __ELECPHYSOPERATORS2D_HPP__
2 #define __ELECPHYSOPERATORS2D_HPP__
3 
4 #include <stdlib.h>
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:
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__
boost::shared_ptr< PreviousData > commonDatav
boost::shared_ptr< PreviousData > commonData
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
structure for User Loop Methods on finite elementsIt can be used to calculate stiffness matrices,...
static FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:141
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
MoFEMErrorCode preProcess()
function is run at the beginning of loop
OpAssembleLhsTauV(std::string flux_field, std::string mass_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
OpSolveRecovery(const std::string &mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Method runge_kutta4)
OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
boost::shared_ptr< PreviousData > commonData
OpEssentialBC(const std::string &flux_field, Range &essential_bd_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::function< double(const double, const double)> FUval
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
boost::shared_ptr< PreviousData > & dataU
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
const double essen_value
boost::shared_ptr< PreviousData > commonData
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
boost::shared_ptr< PreviousData > dataU
boost::function< double(const double, const double)> FUval
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
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)
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
OpAssembleSlowRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &common_datau, boost::shared_ptr< PreviousData > &common_datav, FUVal rhs_u)
boost::shared_ptr< PreviousData > commonDatau
PetscReal ts_t
time
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
boost::function< double(const double, const double)> FUVal
const VectorInt & getIndices() const
Get global indices of dofs on entity.
bool sYmm
If true assume that matrix is symmetric structure.
boost::shared_ptr< PreviousData > commonDatau
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
PetscReal ts_a
shift for U_tt (see PETSc Time Solver)
ForcesAndSourcesCore::UserDataOperator UserDataOperator
MoFEM::FaceElementForcesAndSourcesCore FaceEle
PetscInt ts_step
time step
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
DataForcesAndSourcesCore::EntData EntData
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
auto getFTensor0IntegrationWeight()
Get integration weights.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
MoFEMErrorCode MatSetValues(Mat M, const DataForcesAndSourcesCore::EntData &row_data, const DataForcesAndSourcesCore::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProc > &post_proc)
BoundaryEle::UserDataOperator OpBoundaryEle
MoFEMErrorCode VecSetValues(Vec V, const DataForcesAndSourcesCore::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
#define CHKERR
Inline error check.
Definition: definitions.h:602
boost::shared_ptr< PreviousData > & dataV
OpAssembleStiffRhsTau(std::string flux_field, boost::shared_ptr< PreviousData > &data, std::map< int, BlockData > &block_map)
boost::function< double(const double, const double, const double)> Method
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
boost::shared_ptr< PostProc > postProc
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
Operator for bi-linear form, usually to calculate values on left hand side.
Postprocess on face.
const VectorDofs & getFieldDofs() const
get dofs data stature FEDofEntity
FTensor::Tensor2< FTensor::PackPtr< double *, Tensor_Dim0 *Tensor_Dim1 >, Tensor_Dim0, Tensor_Dim1 > getFTensor2DiffN(FieldApproximationBase base)
Get derivatives of base functions for Hdiv space.
OpAssembleLhsVV(std::string mass_field, boost::shared_ptr< PreviousData > &datau, boost::shared_ptr< PreviousData > &datav, FUval Drhs_u)
boost::shared_ptr< PreviousData > dataV
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Data on single entity (This is passed as argument to DataOperator::doWork)
OpInitialMass(const std::string &mass_field, Range &inner_surface, double &init_val)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorDouble & getFieldData() const
get dofs values
boost::shared_ptr< PreviousData > commonDatav
FaceEle::UserDataOperator OpFaceEle
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FTensor::Index< 'i', 3 > i
OpAssembleLhsTauTau(std::string flux_field, boost::shared_ptr< PreviousData > &commonData, std::map< int, BlockData > &block_map)