v0.9.1
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, 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:
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__
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
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
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:75
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
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.
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
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
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:413
ublas::vector< double, DoubleAllocator > VectorDouble
Definition: Types.hpp:73
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)