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>
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,
149  Range &inner_surface,
150  double init_val)
151  : OpVolEle(mass_field, OpVolEle::OPROW)
152  , innerSurface(inner_surface),
153  initVal(init_val) {}
156  Range &innerSurface;
157  double initVal;
158  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
160  int nb_dofs = data.getFieldData().size();
161  if (nb_dofs) {
162  EntityHandle fe_ent = getFEEntityHandle();
163  bool is_inner_side = (innerSurface.find(fe_ent) != innerSurface.end());
164  if (is_inner_side) {
165  int nb_gauss_pts = getGaussPts().size2();
166  if (nb_dofs != static_cast<int>(data.getN().size2()))
167  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
168  "wrong number of dofs");
169  nN.resize(nb_dofs, nb_dofs, false);
170  nF.resize(nb_dofs, false);
171  nN.clear();
172  nF.clear();
173 
174  auto t_row_mass = data.getFTensor0N();
175  auto t_w = getFTensor0IntegrationWeight();
176  const double vol = getMeasure();
177 
178  for (int gg = 0; gg < nb_gauss_pts; gg++) {
179  const double a = t_w * vol;
180  // double r = ((double)rand() / (RAND_MAX));
181  for (int rr = 0; rr != nb_dofs; rr++) {
182  auto t_col_mass = data.getFTensor0N(gg, 0);
183  nF[rr] += a * initVal * t_row_mass;
184  for (int cc = 0; cc != nb_dofs; cc++) {
185  nN(rr, cc) += a * t_row_mass * t_col_mass;
186  ++t_col_mass;
187  }
188  ++t_row_mass;
189  }
190  ++t_w;
191  }
192 
194  cholesky_solve(nN, nF, ublas::lower());
195 
196  for (auto &dof : data.getFieldDofs()) {
197  dof->getFieldData() = nF[dof->getEntDofIdx()];
198 
199  // this is only to check
200  // data.getFieldData()[dof->getEntDofIdx()] = nF[dof->getEntDofIdx()];
201  }
202  }
203  }
205  }
206 };
207 
208 
209 struct OpSolveRecovery : public OpVolEle {
210  typedef boost::function<double(const double, const double, const double)>
212  OpSolveRecovery(const std::string &mass_field,
213  boost::shared_ptr<PreviousData> &data_u,
214  boost::shared_ptr<PreviousData> &data_v,
215  Method runge_kutta4)
216  : OpVolEle(mass_field, OpVolEle::OPROW)
217  , dataU(data_u)
218  , dataV(data_v)
219  , rungeKutta4(runge_kutta4)
220  {}
221  boost::shared_ptr<PreviousData> dataU;
222  boost::shared_ptr<PreviousData> dataV;
224 
227  double initVal;
228  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
230  int nb_dofs = data.getFieldData().size();
231  if (nb_dofs) {
232  int nb_gauss_pts = getGaussPts().size2();
233  if (nb_dofs != static_cast<int>(data.getN().size2()))
234  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
235  "wrong number of dofs");
236  nN.resize(nb_dofs, nb_dofs, false);
237  nF.resize(nb_dofs, false);
238  nN.clear();
239  nF.clear();
240 
241  auto t_val_u = getFTensor0FromVec(dataU->mass_values);
242  auto t_val_v = getFTensor0FromVec(dataV->mass_values);
243 
244  double dt;
245  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
246 
247  auto t_row_mass = data.getFTensor0N();
248  auto t_w = getFTensor0IntegrationWeight();
249  const double vol = getMeasure();
250 
251  for (int gg = 0; gg < nb_gauss_pts; gg++) {
252  const double a = t_w * vol;
253  const double vn = rungeKutta4(t_val_u, t_val_v, dt);
254  for (int rr = 0; rr != nb_dofs; rr++) {
255  auto t_col_mass = data.getFTensor0N(gg, 0);
256  nF[rr] += a * vn * t_row_mass;
257  for (int cc = 0; cc != nb_dofs; cc++) {
258  nN(rr, cc) += a * t_row_mass * t_col_mass;
259  ++t_col_mass;
260  }
261  ++t_row_mass;
262  }
263  ++t_w;
264  ++t_val_u;
265  ++t_val_v;
266  }
267 
269  cholesky_solve(nN, nF, ublas::lower());
270 
271  for (auto &dof : data.getFieldDofs()) {
272  dof->getFieldData() = nF[dof->getEntDofIdx()];
273 
274  // this is only to check
275  // data.getFieldData()[dof->getEntDofIdx()] = nF[dof->getEntDofIdx()];
276  }
277 
278  }
280  }
281 };
282 
283 // Assembly of RHS for explicit (slow)
284 // part//**************************************
285 
286 // 2. RHS for explicit part of the mass balance equation
288 {
289  typedef boost::function<double(const double, const double)>
291  OpAssembleSlowRhsV(std::string mass_field,
292  boost::shared_ptr<PreviousData> &data_u,
293  boost::shared_ptr<PreviousData> &data_v,
294  Feval_u rhs_u)
295  : OpVolEle(mass_field, OpVolEle::OPROW)
296  , dataU(data_u)
297  , dataV(data_v)
298  , rhsU(rhs_u) {}
299 
304 
305  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
307 
308  const int nb_dofs = data.getIndices().size();
309  if (nb_dofs) {
310  if (nb_dofs != static_cast<int>(data.getN().size2()))
311  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
312  "wrong number of dofs");
313  vecF.resize(nb_dofs, false);
314  mat.resize(nb_dofs, nb_dofs, false);
315  vecF.clear();
316  mat.clear();
317  const int nb_integration_pts = getGaussPts().size2();
318  auto t_val_u = getFTensor0FromVec(dataU->mass_values);
319  auto t_val_v = getFTensor0FromVec(dataV->mass_values);
320 
321 
322  auto t_row_v_base = data.getFTensor0N();
323  auto t_w = getFTensor0IntegrationWeight();
324  const double vol = getMeasure();
325 
326 
327  for (int gg = 0; gg != nb_integration_pts; ++gg) {
328 
329  double rhs = rhsU(t_val_u, t_val_v);
330  const double a = vol * t_w;
331 
332  for (int rr = 0; rr != nb_dofs; ++rr) {
333  auto t_col_v_base = data.getFTensor0N(gg, 0);
334  vecF[rr] += a * rhs * t_row_v_base;
335  for (int cc = 0; cc != nb_dofs; ++cc) {
336  mat(rr, cc) += a * t_row_v_base * t_col_v_base;
337  ++t_col_v_base;
338  }
339  ++t_row_v_base;
340  }
341  ++t_val_u;
342  ++t_val_v;
343  ++t_w;
344 
345  }
347  cholesky_solve(mat, vecF, ublas::lower());
348 
349  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
350  PETSC_TRUE);
351  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
352  ADD_VALUES);
353  }
355  }
356 
357 private:
358  boost::shared_ptr<PreviousData> dataU;
359  boost::shared_ptr<PreviousData> dataV;
362 
363 };
364 
365 // // 5. RHS contribution of the natural boundary condition
366 // struct OpAssembleNaturalBCRhsTau : OpFaceEle // R_tau_2
367 // {
368 // OpAssembleNaturalBCRhsTau(std::string flux_field, Range &natural_bd_ents)
369 // : OpFaceEle(flux_field, OpFaceEle::OPROW),
370 // natural_bd_ents(natural_bd_ents) {}
371 
372 // MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
373 // MoFEMFunctionBegin;
374 // const int nb_dofs = data.getIndices().size();
375 
376 // if (nb_dofs) {
377 // EntityHandle row_side_ent = data.getFieldDofs()[0]->getEnt();
378 
379 // bool is_natural =
380 // (natural_bd_ents.find(row_side_ent) != natural_bd_ents.end());
381 // if (is_natural) {
382 // // cerr << "In NaturalBCRhsTau..." << endl;
383 // vecF.resize(nb_dofs, false);
384 // vecF.clear();
385 // const int nb_integration_pts = getGaussPts().size2();
386 // auto t_tau_base = data.getFTensor1N<3>();
387 
388 // auto dir = getDirection();
389 // FTensor::Tensor1<double, 3> t_normal(-dir[1], dir[0], dir[2]);
390 
391 // auto t_w = getFTensor0IntegrationWeight();
392 
393 // for (int gg = 0; gg != nb_integration_pts; ++gg) {
394 // const double a = t_w;
395 // for (int rr = 0; rr != nb_dofs; ++rr) {
396 // vecF[rr] += (t_tau_base(i) * t_normal(i) * natu_value) * a;
397 // ++t_tau_base;
398 // }
399 // ++t_w;
400 // }
401 // CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
402 // PETSC_TRUE);
403 // CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
404 // ADD_VALUES);
405 // }
406 // }
407 // MoFEMFunctionReturn(0);
408 // }
409 
410 // private:
411 // VectorDouble vecF;
412 // Range natural_bd_ents;
413 // };
414 
415 // Assembly of RHS for the implicit (stiff) part excluding the essential
416 // boundary //**********************************
417 // 3. Assembly of F_tau excluding the essential boundary condition
418 template <int dim>
419 struct OpAssembleStiffRhsTau : OpVolEle // F_tau_1
420 {
421  OpAssembleStiffRhsTau(std::string flux_field,
422  boost::shared_ptr<PreviousData> &data)
423  : OpVolEle(flux_field, OpVolEle::OPROW), commonData(data) {}
424 
425  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
427 
428  const int nb_dofs = data.getIndices().size();
429  if (nb_dofs) {
430 
431 
432  vecF.resize(nb_dofs, false);
433  vecF.clear();
434 
435  const int nb_integration_pts = getGaussPts().size2();
436  auto t_flux_value = getFTensor1FromMat<3>(commonData->flux_values);
437  auto t_mass_value = getFTensor0FromVec(commonData->mass_values);
438  auto t_tau_base = data.getFTensor1N<3>();
439 
440  auto t_tau_grad = data.getFTensor2DiffN<3, 3>();
441 
442  auto t_w = getFTensor0IntegrationWeight();
443  const double vol = getMeasure();
444 
445  for (int gg = 0; gg < nb_integration_pts; ++gg) {
446 
447  const double K_inv = 1. / D_tilde;
448  const double a = vol * t_w;
449  for (int rr = 0; rr < nb_dofs; ++rr) {
450  double div_base =
451  t_tau_grad(0, 0) + t_tau_grad(1, 1) + t_tau_grad(2, 2);
452  vecF[rr] += (K_inv * t_tau_base(i) * t_flux_value(i) -
453  div_base * t_mass_value) *
454  a;
455  ++t_tau_base;
456  ++t_tau_grad;
457  }
458  ++t_flux_value;
459  ++t_mass_value;
460  ++t_w;
461  }
462  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
463  PETSC_TRUE);
464  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
465  ADD_VALUES);
466  }
468  }
469 
470 private:
471  boost::shared_ptr<PreviousData> commonData;
473 };
474 // 4. Assembly of F_v
475 template <int dim>
477 {
478  OpAssembleStiffRhsV(std::string mass_field,
479  boost::shared_ptr<PreviousData> &data_u,
480  Range &stim_region)
481  : OpVolEle(mass_field, OpVolEle::OPROW)
482  , dataU(data_u)
483  , stimRegion(stim_region) {}
484 
485  Range &stimRegion;
486 
487  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
489  const int nb_dofs = data.getIndices().size();
490  // cerr << "In StiffRhsV ..." << endl;
491  if (nb_dofs) {
492 
493  vecF.resize(nb_dofs, false);
494  vecF.clear();
495  const int nb_integration_pts = getGaussPts().size2();
496  auto t_dot_u = getFTensor0FromVec(dataU->mass_dots);
497 
498 
499  auto t_div_u = getFTensor0FromVec(dataU->flux_divs);
500  auto t_row_v_base = data.getFTensor0N();
501  auto t_w = getFTensor0IntegrationWeight();
502  const double vol = getMeasure();
503  const double c_time = getFEMethod()->ts_t;
504 
505  double dt;
506  CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
507 
508  double stim = 0.0;
509 
510  double T = 627.0;
511  double duration = 5.0;
512 
513  if (T - dt < c_time && c_time <= T + duration) {
514  EntityHandle stim_ent = getFEEntityHandle();
515  if (stimRegion.find(stim_ent) != stimRegion.end()) {
516  stim = 40.0;
517  } else {
518  stim = 0.0;
519  }
520  }
521 
522  for (int gg = 0; gg < nb_integration_pts; ++gg) {
523  const double a = vol * t_w;
524  for (int rr = 0; rr < nb_dofs; ++rr) {
525  vecF[rr] += t_row_v_base * (t_dot_u + t_div_u - stim) * a;
526  ++t_row_v_base;
527  }
528  ++t_dot_u;
529  ++t_div_u;
530  ++t_w;
531  }
532  CHKERR VecSetOption(getFEMethod()->ts_F, VEC_IGNORE_NEGATIVE_INDICES,
533  PETSC_TRUE);
534  CHKERR VecSetValues(getFEMethod()->ts_F, data, &*vecF.begin(),
535  ADD_VALUES);
536  }
538  }
539 
540 private:
541  boost::shared_ptr<PreviousData> dataU;
542  boost::shared_ptr<PreviousData> dataV;
544 };
545 
546 // Tangent operator
547 // //**********************************************
548 // 7. Tangent assembly for F_tautau excluding the essential boundary condition
549 template <int dim>
550 struct OpAssembleLhsTauTau : OpVolEle // A_TauTau_1
551 {
552  OpAssembleLhsTauTau(std::string flux_field)
553  : OpVolEle(flux_field, flux_field, OpVolEle::OPROWCOL) {
554  sYmm = true;
555  }
556 
557  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
558  EntityType col_type, EntData &row_data,
559  EntData &col_data) {
561  const int nb_row_dofs = row_data.getIndices().size();
562  const int nb_col_dofs = col_data.getIndices().size();
563 
564  if (nb_row_dofs && nb_col_dofs) {
565 
566  mat.resize(nb_row_dofs, nb_col_dofs, false);
567  mat.clear();
568  const int nb_integration_pts = getGaussPts().size2();
569 
570  auto t_row_tau_base = row_data.getFTensor1N<3>();
571 
572  auto t_w = getFTensor0IntegrationWeight();
573  const double vol = getMeasure();
574 
575  for (int gg = 0; gg != nb_integration_pts; ++gg) {
576  const double a = vol * t_w;
577  const double K_inv = 1. / D_tilde;
578  for (int rr = 0; rr != nb_row_dofs; ++rr) {
579  auto t_col_tau_base = col_data.getFTensor1N<3>(gg, 0);
580  for (int cc = 0; cc != nb_col_dofs; ++cc) {
581  mat(rr, cc) += (K_inv * t_row_tau_base(i) * t_col_tau_base(i)) * a;
582  ++t_col_tau_base;
583  }
584  ++t_row_tau_base;
585  }
586  ++t_w;
587  }
588  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
589  ADD_VALUES);
590  if (row_side != col_side || row_type != col_type) {
591  transMat.resize(nb_col_dofs, nb_row_dofs, false);
592  noalias(transMat) = trans(mat);
593  CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
594  &transMat(0, 0), ADD_VALUES);
595  }
596  }
598  }
599 
600 private:
602 };
603 
604 // 9. Assembly of tangent for F_tau_v excluding the essential bc
605 template <int dim>
606 struct OpAssembleLhsTauV : OpVolEle // E_TauV
607 {
608  OpAssembleLhsTauV(std::string flux_field, std::string mass_field)
609  : OpVolEle(flux_field, mass_field, OpVolEle::OPROWCOL) {
610  sYmm = false;
611  }
612 
613  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
614  EntityType col_type, EntData &row_data,
615  EntData &col_data) {
617  const int nb_row_dofs = row_data.getIndices().size();
618  const int nb_col_dofs = col_data.getIndices().size();
619 
620  if (nb_row_dofs && nb_col_dofs) {
621 
622  mat.resize(nb_row_dofs, nb_col_dofs, false);
623  mat.clear();
624  const int nb_integration_pts = getGaussPts().size2();
625  auto t_w = getFTensor0IntegrationWeight();
626  auto t_row_tau_base = row_data.getFTensor1N<3>();
627 
628  auto t_row_tau_grad = row_data.getFTensor2DiffN<3, 3>();
629  const double vol = getMeasure();
630 
631  for (int gg = 0; gg != nb_integration_pts; ++gg) {
632  const double a = vol * t_w;
633 
634  for (int rr = 0; rr != nb_row_dofs; ++rr) {
635  auto t_col_v_base = col_data.getFTensor0N(gg, 0);
636  for (int cc = 0; cc != nb_col_dofs; ++cc) {
637  double div_row_base = t_row_tau_grad(0, 0) + t_row_tau_grad(1, 1) +
638  t_row_tau_grad(2, 2);
639  mat(rr, cc) += -(div_row_base * t_col_v_base) * a;
640  ++t_col_v_base;
641  }
642  ++t_row_tau_base;
643  ++t_row_tau_grad;
644  }
645  ++t_w;
646  }
647  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
648  ADD_VALUES);
649  }
651  }
652 
653 private:
655 };
656 
657 // 10. Assembly of tangent for F_v_tau
658 struct OpAssembleLhsVTau : OpVolEle // C_VTau
659 {
660  OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
661  : OpVolEle(mass_field, flux_field, OpVolEle::OPROWCOL) {
662  sYmm = false;
663  }
664 
665  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
666  EntityType col_type, EntData &row_data,
667  EntData &col_data) {
669  const int nb_row_dofs = row_data.getIndices().size();
670  const int nb_col_dofs = col_data.getIndices().size();
671 
672  if (nb_row_dofs && nb_col_dofs) {
673  mat.resize(nb_row_dofs, nb_col_dofs, false);
674  mat.clear();
675  const int nb_integration_pts = getGaussPts().size2();
676  auto t_w = getFTensor0IntegrationWeight();
677  auto t_row_v_base = row_data.getFTensor0N();
678  const double vol = getMeasure();
679 
680  for (int gg = 0; gg != nb_integration_pts; ++gg) {
681  const double a = vol * t_w;
682  for (int rr = 0; rr != nb_row_dofs; ++rr) {
683  auto t_col_tau_grad = col_data.getFTensor2DiffN<3, 3>(gg, 0);
684  for (int cc = 0; cc != nb_col_dofs; ++cc) {
685  double div_col_base = t_col_tau_grad(0, 0) + t_col_tau_grad(1, 1) + t_col_tau_grad(2, 2);
686  mat(rr, cc) += (t_row_v_base * div_col_base) * a;
687  ++t_col_tau_grad;
688  }
689  ++t_row_v_base;
690  }
691  ++t_w;
692  }
693  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
694  ADD_VALUES);
695  }
697  }
698 
699 private:
701 };
702 
703 // 11. Assembly of tangent for F_v_v
705 {
706  OpAssembleLhsVV(std::string mass_field)
707  : OpVolEle(mass_field, mass_field, OpVolEle::OPROWCOL) {
708  sYmm = true;
709  }
710 
711  MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type,
712  EntityType col_type, EntData &row_data,
713  EntData &col_data) {
715 
716  const int nb_row_dofs = row_data.getIndices().size();
717  const int nb_col_dofs = col_data.getIndices().size();
718  if (nb_row_dofs && nb_col_dofs) {
719 
720  mat.resize(nb_row_dofs, nb_col_dofs, false);
721  mat.clear();
722  const int nb_integration_pts = getGaussPts().size2();
723 
724  auto t_row_v_base = row_data.getFTensor0N();
725 
726  auto t_w = getFTensor0IntegrationWeight();
727  const double ts_a = getFEMethod()->ts_a;
728  const double vol = getMeasure();
729 
730  for (int gg = 0; gg != nb_integration_pts; ++gg) {
731  const double a = vol * t_w;
732 
733  for (int rr = 0; rr != nb_row_dofs; ++rr) {
734  auto t_col_v_base = col_data.getFTensor0N(gg, 0);
735  for (int cc = 0; cc != nb_col_dofs; ++cc) {
736  mat(rr, cc) += (ts_a * t_row_v_base * t_col_v_base) * a;
737 
738  ++t_col_v_base;
739  }
740  ++t_row_v_base;
741  }
742  ++t_w;
743  }
744  CHKERR MatSetValues(getFEMethod()->ts_B, row_data, col_data, &mat(0, 0),
745  ADD_VALUES);
746  if (row_side != col_side || row_type != col_type) {
747  transMat.resize(nb_col_dofs, nb_row_dofs, false);
748  noalias(transMat) = trans(mat);
749  CHKERR MatSetValues(getFEMethod()->ts_B, col_data, row_data,
750  &transMat(0, 0), ADD_VALUES);
751  }
752  }
754  }
755 
756 private:
758 };
759 
760 // struct OpError : public OpFaceEle {
761 // typedef boost::function<double(const double, const double, const double)>
762 // FVal;
763 // typedef boost::function<FTensor::Tensor1<double, 3>(
764 // const double, const double, const double)>
765 // FGrad;
766 // double &eRror;
767 // OpError(FVal exact_value, FVal exact_lap, FGrad exact_grad,
768 // boost::shared_ptr<PreviousData> &prev_data,
769 // std::map<int, BlockData> &block_map, double &err)
770 // : OpFaceEle("ERROR", OpFaceEle::OPROW), exactVal(exact_value),
771 // exactLap(exact_lap), exactGrad(exact_grad), prevData(prev_data),
772 // setOfBlock(block_map), eRror(err) {}
773 // MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
774 // MoFEMFunctionBegin;
775 // const int nb_dofs = data.getFieldData().size();
776 // // cout << "nb_error_dofs : " << nb_dofs << endl;
777 // if (nb_dofs) {
778 // auto find_block_data = [&]() {
779 // EntityHandle fe_ent = getFEEntityHandle();
780 // BlockData *block_raw_ptr = nullptr;
781 // for (auto &m : setOfBlock) {
782 // if (m.second.block_ents.find(fe_ent) != m.second.block_ents.end()) {
783 // block_raw_ptr = &m.second;
784 // break;
785 // }
786 // }
787 // return block_raw_ptr;
788 // };
789 
790 // auto block_data_ptr = find_block_data();
791 // if (!block_data_ptr)
792 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Block not found");
793 
794 // auto &block_data = *block_data_ptr;
795 
796 // auto t_flux_value = getFTensor1FromMat<3>(prevData->flux_values);
797 // // auto t_mass_dot = getFTensor0FromVec(prevData->mass_dots);
798 // auto t_mass_value = getFTensor0FromVec(prevData->mass_values);
799 // // auto t_flux_div = getFTensor0FromVec(prevData->flux_divs);
800 // data.getFieldData().clear();
801 // const double vol = getMeasure();
802 // const int nb_integration_pts = getGaussPts().size2();
803 // auto t_w = getFTensor0IntegrationWeight();
804 // double dt;
805 // CHKERR TSGetTimeStep(getFEMethod()->ts, &dt);
806 // double ct = getFEMethod()->ts_t - dt;
807 // auto t_coords = getFTensor1CoordsAtGaussPts();
808 
809 // FTensor::Tensor1<double, 3> t_exact_flux, t_flux_error;
810 
811 // for (int gg = 0; gg != nb_integration_pts; ++gg) {
812 // const double a = vol * t_w;
813 // double mass_exact = exactVal(t_coords(NX), t_coords(NY), ct);
814 // // double flux_lap = - block_data.B0 * exactLap(t_coords(NX),
815 // // t_coords(NY), ct);
816 // t_exact_flux(i) =
817 // -block_data.B0 * exactGrad(t_coords(NX), t_coords(NY), ct)(i);
818 // t_flux_error(0) = t_flux_value(0) - t_exact_flux(0);
819 // t_flux_error(1) = t_flux_value(1) - t_exact_flux(1);
820 // t_flux_error(2) = t_flux_value(2) - t_exact_flux(2);
821 // double local_error = pow(mass_exact - t_mass_value, 2) +
822 // t_flux_error(i) * t_flux_error(i);
823 // // cout << "flux_div : " << t_flux_div << " flux_exact : " <<
824 // // flux_exact << endl;
825 // data.getFieldData()[0] += a * local_error;
826 // eRror += a * local_error;
827 
828 // ++t_w;
829 // ++t_mass_value;
830 // // ++t_flux_div;
831 // ++t_flux_value;
832 // // ++t_mass_dot;
833 // ++t_coords;
834 // }
835 
836 // data.getFieldDofs()[0]->getFieldData() = data.getFieldData()[0];
837 // }
838 // MoFEMFunctionReturn(0);
839 // }
840 
841 // private:
842 // FVal exactVal;
843 // FVal exactLap;
844 // FGrad exactGrad;
845 // boost::shared_ptr<PreviousData> prevData;
846 // std::map<int, BlockData> setOfBlock;
847 
848 // FTensor::Number<0> NX;
849 // FTensor::Number<1> NY;
850 // };
851 
852 
853 
854 struct Monitor : public FEMethod {
855  Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj<DM> &dm,
856  boost::shared_ptr<PostProcVolumeOnRefinedMesh> &post_proc)
857  : cOmm(comm), rAnk(rank), dM(dm), postProc(post_proc){};
858  MoFEMErrorCode preProcess() { return 0; }
859  MoFEMErrorCode operator()() { return 0; }
862  if (ts_step % save_every_nth_step == 0) {
864  CHKERR postProc->writeFile(
865  "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
866  }
867 
869  }
870 
871 private:
873 
874  boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProc;
875  MPI_Comm cOmm;
876  const int rAnk;
877 };
878 
879 }; // namespace ReactionDiffusion
880 
881 #endif //__ELECPHYSOPERATORS_HPP__
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:143
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
OpSolveRecovery(const std::string &mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Method runge_kutta4)
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
boost::shared_ptr< PreviousData > dataV
auto post_proc
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)> Feval_u
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.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
const double essen_value
OpInitialMass(const std::string &mass_field, Range &inner_surface, double init_val)
boost::shared_ptr< PreviousData > dataU
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
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
const MatrixAdaptor getVectorN(const FieldApproximationBase base, const int gg)
get Hdiv of base functions at Gauss pts
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
PetscReal ts_t
time
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
const VectorInt & getIndices() const
Get global indices of dofs on entity.
OpAssembleLhsTauV(std::string flux_field, std::string mass_field)
OpAssembleLhsVV(std::string mass_field)
bool sYmm
If true assume that matrix is symmetric structure.
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
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
MoFEM::FaceElementForcesAndSourcesCore FaceEle
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
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.
OpAssembleSlowRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Feval_u rhs_u)
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
#define CHKERR
Inline error check.
Definition: definitions.h:604
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcVolumeOnRefinedMesh > &post_proc)
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.
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.
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.
OpAssembleStiffRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &data_u, Range &stim_region)
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)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:74
boost::shared_ptr< PreviousData > dataV
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
boost::shared_ptr< PreviousData > dataU
const VectorDouble & getFieldData() const
get dofs values
DataForcesAndSourcesCore::EntData EntData
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1N(FieldApproximationBase base)
Get base functions for Hdiv/Hcurl spaces.
OpEssentialBC(std::string flux_field, Range &essential_bd_ents)
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
FTensor::Index< 'i', 3 > i
OpAssembleStiffRhsTau(std::string flux_field, boost::shared_ptr< PreviousData > &data)