v0.13.1
ElecPhysOperators.hpp
Go to the documentation of this file.
1#ifndef __ELECPHYSOPERATORS_HPP__
2#define __ELECPHYSOPERATORS_HPP__
3
4#include <stdlib.h>
5#include <BasicFiniteElements.hpp>
6
8
10
12
15
17
18
19
20const int save_every_nth_step = 1;
21
22const double essen_value = 0;
23
25
26
27// problem parameters
28const double alpha = 0.01;
29const double gma = 0.002;
30const double b = 0.15;
31const double c = 8.00;
32const double mu1 = 0.20;
33const double mu2 = 0.30;
34
35const double D_tilde = 0.01;
36
37
38
42
45
47
50
52 jac.resize(3, 3, false);
53 inv_jac.resize(3, 3, false);
54 }
55};
56
57
58struct OpEssentialBC : public OpFaceEle {
59 OpEssentialBC(std::string flux_field, Range &essential_bd_ents)
61 }
64 int nb_dofs = data.getIndices().size();
65
66 if (nb_dofs) {
67 EntityHandle fe_ent = getFEEntityHandle();
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
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
135private:
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
147struct 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) {}
155 double initVal;
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
207struct 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;
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
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
355private:
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
416template <int dim>
418{
419 OpAssembleStiffRhsTau(std::string flux_field,
420 boost::shared_ptr<PreviousData> &data)
421 : OpVolEle(flux_field, OpVolEle::OPROW), commonData(data) {}
422
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
468private:
469 boost::shared_ptr<PreviousData> commonData;
471};
472// 4. Assembly of F_v
473template <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
484
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
538private:
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
547template <int dim>
548struct 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
598private:
600};
601
602// 9. Assembly of tangent for F_tau_v excluding the essential bc
603template <int dim>
604struct 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
651private:
653};
654
655// 10. Assembly of tangent for F_v_tau
656struct 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
697private:
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
754private:
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
852struct 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){};
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
869private:
870 SmartPetscObj<DM> dM;
871
872 boost::shared_ptr<PostProcVolumeOnRefinedMesh> postProc;
873 MPI_Comm cOmm;
874 const int rAnk;
875};
876
877}; // namespace ReactionDiffusion
878
879#endif //__ELECPHYSOPERATORS_HPP__
constexpr double a
void cholesky_solve(const TRIA &L, VEC &x, ublas::lower)
solve system L L^T x = b inplace
Definition: cholesky.hpp:221
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
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ HVEC0
Definition: definitions.h:199
@ HVEC1
Definition: definitions.h:199
@ HVEC2
Definition: definitions.h:199
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:545
double dt
Definition: heat_method.cpp:40
const double T
DataForcesAndSourcesCore::EntData EntData
const double essen_value
FTensor::Index< 'i', 3 > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
UBlasVector< double > VectorDouble
Definition: Types.hpp:79
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
double duration
impulse duration (ns)
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcVolumeOnRefinedMesh > &post_proc)
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsTauV(std::string flux_field, std::string mass_field)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsVTau(std::string mass_field, std::string flux_field)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, int col_side, EntityType row_type, EntityType col_type, EntData &row_data, EntData &col_data)
OpAssembleLhsVV(std::string mass_field)
OpAssembleSlowRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Feval_u rhs_u)
boost::function< double(const double, const double)> Feval_u
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > dataU
boost::shared_ptr< PreviousData > dataV
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::shared_ptr< PreviousData > commonData
OpAssembleStiffRhsTau(std::string flux_field, boost::shared_ptr< PreviousData > &data)
boost::shared_ptr< PreviousData > dataV
boost::shared_ptr< PreviousData > dataU
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpAssembleStiffRhsV(std::string mass_field, boost::shared_ptr< PreviousData > &data_u, Range &stim_region)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpEssentialBC(std::string flux_field, Range &essential_bd_ents)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
OpInitialMass(const std::string &mass_field, Range &inner_surface, double init_val)
OpSolveRecovery(const std::string &mass_field, boost::shared_ptr< PreviousData > &data_u, boost::shared_ptr< PreviousData > &data_v, Method runge_kutta4)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
boost::function< double(const double, const double, const double)> Method
boost::shared_ptr< PreviousData > dataV
boost::shared_ptr< PreviousData > dataU
bool sYmm
If true assume that matrix is symmetric structure.
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
@ OPROWCOL
operator doWork is executed on FE rows &columns
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
PetscReal ts_t
time
PetscReal ts_a
shift for U_t (see PETSc Time Solver)