v0.14.0
Loading...
Searching...
No Matches
ElecPhysOperators.hpp
Go to the documentation of this file.
1#ifndef __ELECPHYSOPERATORS_HPP__
2#define __ELECPHYSOPERATORS_HPP__
3
4#include <stdlib.h>
6
8
10
12
15
16using EntData = DataForcesAndSourcesCore::EntData;
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
40 MatrixDouble flux_values;
41 VectorDouble flux_divs;
42
43 VectorDouble mass_dots;
44 VectorDouble mass_values;
45
46 VectorDouble slow_values;
47
48 MatrixDouble jac;
49 MatrixDouble inv_jac;
50
52 jac.resize(3, 3, false);
53 inv_jac.resize(3, 3, false);
54 }
55};
56
57
58struct OpEssentialBC : public OpFaceEle {
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
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:
136 MatrixDouble nN;
137 VectorDouble nF;
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) {}
152 MatrixDouble nN;
153 VectorDouble nF;
155 double initVal;
156 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
158 int nb_dofs = data.getFieldData().size();
159 if (nb_dofs) {
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
223 MatrixDouble nN;
224 VectorDouble nF;
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
355private:
356 boost::shared_ptr<PreviousData> dataU;
357 boost::shared_ptr<PreviousData> dataV;
358 VectorDouble vecF;
359 MatrixDouble mat;
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
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
468private:
469 boost::shared_ptr<PreviousData> commonData;
470 VectorDouble vecF;
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
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
538private:
539 boost::shared_ptr<PreviousData> dataU;
540 boost::shared_ptr<PreviousData> dataV;
541 VectorDouble vecF;
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:
599 MatrixDouble mat, transMat;
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:
652 MatrixDouble mat;
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:
698 MatrixDouble mat;
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:
755 MatrixDouble mat, transMat;
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:
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 ...
@ HVEC0
@ HVEC1
@ HVEC2
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:572
double dt
DataForcesAndSourcesCore::EntData EntData
FTensor::Index< 'i', 3 > i
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
MoFEMErrorCode MatSetValues(Mat M, const EntitiesFieldData::EntData &row_data, const EntitiesFieldData::EntData &col_data, const double *ptr, InsertMode iora)
Assemble PETSc matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Monitor(MPI_Comm &comm, const int &rank, SmartPetscObj< DM > &dm, boost::shared_ptr< PostProcVolumeOnRefinedMesh > &post_proc)
boost::shared_ptr< PostProcVolumeOnRefinedMesh > postProc
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode postProcess()
function is run at the end of loop
MoFEMErrorCode preProcess()
function is run at the beginning of loop
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)
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.
structure for User Loop Methods on finite elements
MatrixDouble & getNormalsAtGaussPts()
if higher order geometry return normals at Gauss pts.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ 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
intrusive_ptr for managing petsc objects
PetscReal ts_t
time
PetscReal ts_a
shift for U_t (see PETSc Time Solver)
PetscInt ts_step
time step number