v0.13.1
Loading...
Searching...
No Matches
FreeSurfaceOps.hpp
Go to the documentation of this file.
1namespace FreeSurfaceOps {
2
4 OpCalculateLift(const std::string field_name,
5 boost::shared_ptr<VectorDouble> p_ptr,
6 boost::shared_ptr<VectorDouble> lift_ptr,
7 boost::shared_ptr<Range> ents_ptr)
9 pPtr(p_ptr), liftPtr(lift_ptr), entsPtr(ents_ptr) {
10 std::fill(&doEntities[MBVERTEX], &doEntities[MBMAXTYPE], false);
11 doEntities[MBEDGE] = true;
12 }
13
14 MoFEMErrorCode doWork(int row_side, EntityType row_type,
15 HookeElement::EntData &data) {
17
18 const auto fe_ent = getFEEntityHandle();
19 if (entsPtr->find(fe_ent) != entsPtr->end()) {
20
22 auto t_p = getFTensor0FromVec(*pPtr);
23 auto t_normal = getFTensor1Normal();
24 auto t_coords = getFTensor1CoordsAtGaussPts();
25 auto t_lift = getFTensor1FromArray<SPACE_DIM, SPACE_DIM>(*liftPtr);
26
27 const auto nb_int_points = getGaussPts().size2();
28
29 for (int gg = 0; gg != nb_int_points; gg++) {
30
31 const double r = t_coords(0);
32 const double alpha = cylindrical(r) * t_w / 2;
33 t_lift(i) -= t_normal(i) * (t_p * alpha);
34
35 ++t_w;
36 ++t_p;
37 ++t_coords;
38 }
39 }
40
42 }
43
44private:
45 boost::shared_ptr<VectorDouble> pPtr;
46 boost::shared_ptr<VectorDouble> liftPtr;
47 boost::shared_ptr<Range> entsPtr;
48};
49
52 boost::shared_ptr<MatrixDouble> u_ptr)
55 uPtr(u_ptr) {}
56
57 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
59
60 auto t_w = getFTensor0IntegrationWeight();
61 auto t_normal = getFTensor1Normal();
62 auto t_u = getFTensor1FromMat<SPACE_DIM>(*uPtr);
63 auto t_row_base = row_data.getFTensor0N();
64 auto t_coords = getFTensor1CoordsAtGaussPts();
65
66 for (int gg = 0; gg != nbIntegrationPts; gg++) {
67
68 const double r = t_coords(0);
69 const double alpha = t_w * cylindrical(r);
70
71 int bb = 0;
72 for (; bb != nbRows; ++bb) {
73 locF[bb] += alpha * t_row_base * (t_normal(i) * t_u(i));
74 ++t_row_base;
75 }
76
77 for (; bb < nbRowBaseFunctions; ++bb)
78 ++t_row_base;
79
80 ++t_w;
81 ++t_u;
82 ++t_coords;
83 }
84
86 }
87
88private:
89 boost::shared_ptr<MatrixDouble> uPtr;
90};
91
93 OpNormalForceRhs(const std::string field_name,
94 boost::shared_ptr<VectorDouble> lambda_ptr)
96 AssemblyDomainEleOp::OPROW),
97 lambdaPtr(lambda_ptr) {}
98
99 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data) {
101
102 auto t_w = getFTensor0IntegrationWeight();
103 auto t_normal = getFTensor1Normal();
104 auto t_lambda = getFTensor0FromVec(*lambdaPtr);
105 auto t_row_base = row_data.getFTensor0N();
106 auto t_coords = getFTensor1CoordsAtGaussPts();
107
108 for (int gg = 0; gg != nbIntegrationPts; gg++) {
109
110 auto t_nf = getFTensor1FromArray<U_FIELD_DIM, U_FIELD_DIM>(locF);
111
112 const double r = t_coords(0);
113 const double alpha = t_w * cylindrical(r);
114
115 int bb = 0;
116 for (; bb != nbRows / U_FIELD_DIM; ++bb) {
117
118 t_nf(i) += alpha * t_row_base * t_normal(i) * t_lambda;
119 ++t_row_base;
120 ++t_nf;
121 }
122
123 for (; bb < nbRowBaseFunctions; ++bb)
124 ++t_row_base;
125
126 ++t_w;
127 ++t_lambda;
128 ++t_coords;
129 }
130
132 }
133
134private:
135 boost::shared_ptr<VectorDouble> lambdaPtr;
136};
137
139
140 OpNormalConstrainLhs(const std::string field_name_row,
141 const std::string field_name_col)
142 : AssemblyBoundaryEleOp(field_name_row, field_name_col,
143 AssemblyBoundaryEleOp::OPROWCOL) {
144 assembleTranspose = true;
145 sYmm = false;
146 }
147
148 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
149 EntitiesFieldData::EntData &col_data) {
151
152 auto t_w = getFTensor0IntegrationWeight();
153 auto t_normal = getFTensor1Normal();
154 auto t_row_base = row_data.getFTensor0N();
155 auto t_coords = getFTensor1CoordsAtGaussPts();
156
157 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
158
159 auto t_mat = getFTensor1FromPtr<U_FIELD_DIM>(&locMat(0, 0));
160
161 const double r = t_coords(0);
162 const double alpha = t_w * cylindrical(r);
163
164 int rr = 0;
165 for (; rr != nbRows; ++rr) {
166
167 auto t_col_base = col_data.getFTensor0N(gg, 0);
168 const double a = alpha * t_row_base;
169
170 for (int cc = 0; cc != nbCols / U_FIELD_DIM; ++cc) {
171 t_mat(i) += (a * t_col_base) * t_normal(i);
172 ++t_col_base;
173 ++t_mat;
174 }
175 ++t_row_base;
176 }
177
178 for (; rr < nbRowBaseFunctions; ++rr)
179 ++t_row_base;
180
181 ++t_w;
182 ++t_coords;
183 }
184
186 };
187};
188
189/**
190 * @brief Rhs for U
191 *
192 */
194
195 OpRhsU(const std::string field_name,
196 boost::shared_ptr<MatrixDouble> dot_u_ptr,
197 boost::shared_ptr<MatrixDouble> u_ptr,
198 boost::shared_ptr<MatrixDouble> grad_u_ptr,
199 boost::shared_ptr<VectorDouble> h_ptr,
200 boost::shared_ptr<MatrixDouble> grad_h_ptr,
201 boost::shared_ptr<VectorDouble> g_ptr,
202 boost::shared_ptr<VectorDouble> p_ptr)
204 dotUPtr(dot_u_ptr), uPtr(u_ptr), gradUPtr(grad_u_ptr), hPtr(h_ptr),
205 gradHPtr(grad_h_ptr), gPtr(g_ptr), pPtr(p_ptr) {}
206
207 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
209
210 const double vol = getMeasure();
211 auto t_dot_u = getFTensor1FromMat<U_FIELD_DIM>(*dotUPtr);
212 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
213 auto t_p = getFTensor0FromVec(*pPtr);
214 auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*gradUPtr);
215 auto t_h = getFTensor0FromVec(*hPtr);
216 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
217 auto t_g = getFTensor0FromVec(*gPtr);
218 auto t_coords = getFTensor1CoordsAtGaussPts();
219
220 auto t_base = data.getFTensor0N();
221 auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
222
223 auto t_w = getFTensor0IntegrationWeight();
224
232
233 t_buoyancy(i) = 0;
234
235 for (int gg = 0; gg != nbIntegrationPts; gg++) {
236
237 const double r = t_coords(0);
238 const double alpha = t_w * vol * cylindrical(r);
239
240 const double rho = phase_function(t_h, rho_diff, rho_ave);
241 const double mu = phase_function(t_h, mu_diff, mu_ave);
242
243 auto t_D = get_D(2 * mu);
244
245 t_inertia_force(i) = (rho * alpha) * (t_dot_u(i));
246 t_buoyancy(SPACE_DIM - 1) = -(alpha * rho * a0) * t_h;
247 t_phase_force(i) = -alpha * kappa * t_g * t_grad_h(i);
248 t_convection(i) = (rho * alpha) * (t_u(j) * t_grad_u(i, j));
249
250 t_stress(i, j) =
251 alpha * (t_D(i, j, k, l) * t_grad_u(k, l) + t_kd(i, j) * t_p);
252
253 auto t_nf = getFTensor1FromArray<U_FIELD_DIM, U_FIELD_DIM>(locF);
254
255 t_forces(i) = t_inertia_force(i) + t_buoyancy(i) + t_convection(i) +
256 t_phase_force(i);
257
258 int bb = 0;
259 for (; bb != nbRows / U_FIELD_DIM; ++bb) {
260
261 t_nf(i) += t_base * t_forces(i);
262 t_nf(i) += t_diff_base(j) * t_stress(i, j);
263
264 // When we move to C++17 add if constexpr()
265 if constexpr (coord_type == CYLINDRICAL) {
266 t_nf(0) += (t_base * (alpha / t_coords(0))) * (2 * mu * t_u(0) + t_p);
267 }
268
269 ++t_base;
270 ++t_diff_base;
271 ++t_nf;
272 }
273
274 for (; bb < nbRowBaseFunctions; ++bb) {
275 ++t_diff_base;
276 ++t_base;
277 }
278
279 ++t_dot_u;
280 ++t_u;
281 ++t_grad_u;
282 ++t_h;
283 ++t_grad_h;
284 ++t_g;
285 ++t_p;
286
287 ++t_w;
288 ++t_coords;
289 }
290
292 }
293
294private:
295 boost::shared_ptr<MatrixDouble> dotUPtr;
296 boost::shared_ptr<MatrixDouble> uPtr;
297 boost::shared_ptr<MatrixDouble> gradUPtr;
298 boost::shared_ptr<VectorDouble> hPtr;
299 boost::shared_ptr<MatrixDouble> gradHPtr;
300 boost::shared_ptr<VectorDouble> gPtr;
301 boost::shared_ptr<VectorDouble> pPtr;
302};
303
304/**
305 * @brief Lhs for U dU
306 *
307 */
309
310 OpLhsU_dU(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
311 boost::shared_ptr<MatrixDouble> grad_u_ptr,
312 boost::shared_ptr<VectorDouble> h_ptr)
314 AssemblyDomainEleOp::OPROWCOL),
315 uPtr(u_ptr), gradUPtr(grad_u_ptr), hPtr(h_ptr) {
316 sYmm = false;
317 assembleTranspose = false;
318 }
319
320 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
321 EntitiesFieldData::EntData &col_data) {
323
324 const double vol = getMeasure();
325 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
326 auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*gradUPtr);
327 auto t_h = getFTensor0FromVec(*hPtr);
328 auto t_coords = getFTensor1CoordsAtGaussPts();
329
330 auto t_row_base = row_data.getFTensor0N();
331 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
332
333 auto t_w = getFTensor0IntegrationWeight();
334
335 auto get_mat = [&](const int rr) {
336 return getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(locMat, rr);
337 };
338
339 auto ts_a = getTSa();
341
342 for (int gg = 0; gg != nbIntegrationPts; gg++) {
343
344 const double r = t_coords(0);
345 const double alpha = t_w * vol * cylindrical(r);
346 const double rho = phase_function(t_h, rho_diff, rho_ave);
347 const double mu = phase_function(t_h, mu_diff, mu_ave);
348
349 const double beta0 = alpha * rho;
350 const double beta1 = beta0 * ts_a;
351 auto t_D = get_D(2 * mu);
352
353 int rr = 0;
354 for (; rr != nbRows / U_FIELD_DIM; ++rr) {
355
356 auto t_mat = get_mat(rr * U_FIELD_DIM);
357 auto t_col_base = col_data.getFTensor0N(gg, 0);
358 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
359
361 // I mix up the indices here so that it behaves like a
362 // Dg. That way I don't have to have a separate wrapper
363 // class Christof_Expr, which simplifies things.
364 t_d_stress(l, j, k) = t_D(i, j, k, l) * (alpha * t_row_diff_base(i));
365
366 for (int cc = 0; cc != nbCols / U_FIELD_DIM; ++cc) {
367
368 const double bb = t_row_base * t_col_base;
369
370 t_mat(i, j) += (beta1 * bb) * t_kd(i, j);
371 t_mat(i, j) += (beta0 * bb) * t_grad_u(i, j);
372 t_mat(i, j) +=
373 (beta0 * t_row_base) * t_kd(i, j) * (t_col_diff_base(k) * t_u(k));
374 t_mat(i, j) += t_d_stress(i, j, k) * t_col_diff_base(k);
375
376 // When we move to C++17 add if constexpr()
377 if constexpr (coord_type == CYLINDRICAL) {
378 t_mat(0, 0) += (bb * (alpha / t_coords(0))) * (2 * mu);
379 }
380
381 ++t_mat;
382 ++t_col_base;
383 ++t_col_diff_base;
384 }
385
386 ++t_row_base;
387 ++t_row_diff_base;
388 }
389
390 for (; rr < nbRowBaseFunctions; ++rr) {
391 ++t_row_diff_base;
392 ++t_row_base;
393 }
394
395 ++t_u;
396 ++t_grad_u;
397 ++t_h;
398
399 ++t_coords;
400 ++t_w;
401 }
402
404 }
405
406private:
407 boost::shared_ptr<MatrixDouble> uPtr;
408 boost::shared_ptr<MatrixDouble> gradUPtr;
409 boost::shared_ptr<VectorDouble> hPtr;
410};
411
412/**
413 * @brief Lhs for U dH
414 *
415 */
417
418 OpLhsU_dH(const std::string field_name_u, const std::string field_name_h,
419 boost::shared_ptr<MatrixDouble> dot_u_ptr,
420 boost::shared_ptr<MatrixDouble> u_ptr,
421 boost::shared_ptr<MatrixDouble> grad_u_ptr,
422 boost::shared_ptr<VectorDouble> h_ptr,
423 boost::shared_ptr<VectorDouble> g_ptr)
424 : AssemblyDomainEleOp(field_name_u, field_name_h,
425 AssemblyDomainEleOp::OPROWCOL),
426 dotUPtr(dot_u_ptr), uPtr(u_ptr), gradUPtr(grad_u_ptr), hPtr(h_ptr),
427 gPtr(g_ptr) {
428 sYmm = false;
429 assembleTranspose = false;
430 }
431
432 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
433 EntitiesFieldData::EntData &col_data) {
435
436 const double vol = getMeasure();
437 auto t_dot_u = getFTensor1FromMat<U_FIELD_DIM>(*dotUPtr);
438 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
439 auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*gradUPtr);
440 auto t_h = getFTensor0FromVec(*hPtr);
441 auto t_g = getFTensor0FromVec(*gPtr);
442 auto t_coords = getFTensor1CoordsAtGaussPts();
443
444 auto t_row_base = row_data.getFTensor0N();
445 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
446
447 auto t_w = getFTensor0IntegrationWeight();
448
451 FTensor::Tensor1<double, U_FIELD_DIM> t_inertia_force_dh;
455
456 t_buoyancy_dh(i) = 0;
457
458 for (int gg = 0; gg != nbIntegrationPts; gg++) {
459
460 const double r = t_coords(0);
461 const double alpha = t_w * vol * cylindrical(r);
462
463 const double rho = phase_function(t_h, rho_diff, rho_ave);
464 const double rho_dh = d_phase_function_h(t_h, rho_diff);
465 const double mu_dh = d_phase_function_h(t_h, mu_diff);
466
467 auto t_D_dh = get_D(2 * mu_dh);
468
469 t_inertia_force_dh(i) = (alpha * rho_dh) * t_dot_u(i);
470 t_buoyancy_dh(SPACE_DIM - 1) = -(alpha * a0) * (rho + rho_dh * t_h);
471 t_convection_dh(i) = (rho_dh * alpha) * (t_u(j) * t_grad_u(i, j));
472 const double t_phase_force_g_dh = -alpha * kappa * t_g;
473 t_forces_dh(i) =
474 t_inertia_force_dh(i) + t_buoyancy_dh(i) + t_convection_dh(i);
475
476 t_stress_dh(i, j) = alpha * (t_D_dh(i, j, k, l) * t_grad_u(k, l));
477
478 int rr = 0;
479 for (; rr != nbRows / U_FIELD_DIM; ++rr) {
480
481 auto t_mat =
482 getFTensor1FromMat<U_FIELD_DIM, 1>(locMat, rr * U_FIELD_DIM);
483 auto t_col_base = col_data.getFTensor0N(gg, 0);
484 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
485
486 for (int cc = 0; cc != nbCols; ++cc) {
487
488 const double bb = t_row_base * t_col_base;
489 t_mat(i) += t_forces_dh(i) * bb;
490 t_mat(i) += (t_phase_force_g_dh * t_row_base) * t_col_diff_base(i);
491 t_mat(i) += (t_row_diff_base(j) * t_col_base) * t_stress_dh(i, j);
492
493 // When we move to C++17 add if constexpr()
494 if constexpr (coord_type == CYLINDRICAL) {
495 t_mat(0) += (bb * (alpha / t_coords(0))) * (2 * mu_dh * t_u(0));
496 }
497
498 ++t_mat;
499 ++t_col_base;
500 ++t_col_diff_base;
501 }
502
503 ++t_row_base;
504 ++t_row_diff_base;
505 }
506
507 for (; rr < nbRowBaseFunctions; ++rr) {
508 ++t_row_diff_base;
509 ++t_row_base;
510 }
511
512 ++t_dot_u;
513 ++t_u;
514 ++t_grad_u;
515 ++t_h;
516 ++t_g;
517 ++t_coords;
518 ++t_w;
519 }
520
522 }
523
524private:
525 boost::shared_ptr<MatrixDouble> dotUPtr;
526 boost::shared_ptr<MatrixDouble> uPtr;
527 boost::shared_ptr<MatrixDouble> gradUPtr;
528 boost::shared_ptr<VectorDouble> hPtr;
529 boost::shared_ptr<VectorDouble> gPtr;
530};
531
532/**
533 * @brief Lhs for G dH
534 *
535 */
537
538 OpLhsU_dG(const std::string field_name_u, const std::string field_name_h,
539 boost::shared_ptr<MatrixDouble> grad_h_ptr)
540 : AssemblyDomainEleOp(field_name_u, field_name_h,
541 AssemblyDomainEleOp::OPROWCOL),
542 gradHPtr(grad_h_ptr) {
543 sYmm = false;
544 assembleTranspose = false;
545 }
546
547 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
548 EntitiesFieldData::EntData &col_data) {
550
551 const double vol = getMeasure();
552 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
553 auto t_coords = getFTensor1CoordsAtGaussPts();
554
555 auto t_row_base = row_data.getFTensor0N();
556 auto t_w = getFTensor0IntegrationWeight();
557
558 for (int gg = 0; gg != nbIntegrationPts; gg++) {
559
560 const double r = t_coords(0);
561 const double alpha = t_w * vol * cylindrical(r);
562
564 t_phase_force_dg(i) = -alpha * kappa * t_grad_h(i);
565
566 int rr = 0;
567 for (; rr != nbRows / U_FIELD_DIM; ++rr) {
568 auto t_mat =
569 getFTensor1FromMat<U_FIELD_DIM, 1>(locMat, rr * U_FIELD_DIM);
570 auto t_col_base = col_data.getFTensor0N(gg, 0);
571
572 for (int cc = 0; cc != nbCols; ++cc) {
573 const double bb = t_row_base * t_col_base;
574 t_mat(i) += t_phase_force_dg(i) * bb;
575
576 ++t_mat;
577 ++t_col_base;
578 }
579
580 ++t_row_base;
581 }
582
583 for (; rr < nbRowBaseFunctions; ++rr) {
584 ++t_row_base;
585 }
586
587 ++t_grad_h;
588 ++t_coords;
589 ++t_w;
590 }
591
593 }
594
595private:
596 boost::shared_ptr<MatrixDouble> gradHPtr;
597};
598
599template <bool I> struct OpRhsH : public AssemblyDomainEleOp {
600
601 OpRhsH(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
602 boost::shared_ptr<VectorDouble> dot_h_ptr,
603 boost::shared_ptr<VectorDouble> h_ptr,
604 boost::shared_ptr<MatrixDouble> grad_h_ptr,
605 boost::shared_ptr<MatrixDouble> grad_g_ptr)
607 uPtr(u_ptr), dotHPtr(dot_h_ptr), hPtr(h_ptr), gradHPtr(grad_h_ptr),
608 gradGPtr(grad_g_ptr) {}
609
610 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
612
613 const double vol = getMeasure();
614 auto t_w = getFTensor0IntegrationWeight();
615 auto t_coords = getFTensor1CoordsAtGaussPts();
616 auto t_base = data.getFTensor0N();
617 auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
618
619#ifndef NDEBUG
620 if(data.getDiffN().size1() != data.getN().size1())
621 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
622 if (data.getDiffN().size2() != data.getN().size2() * SPACE_DIM) {
623 MOFEM_LOG("SELF", Sev::error)
624 << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
625 MOFEM_LOG("SELF", Sev::error) << data.getN();
626 MOFEM_LOG("SELF", Sev::error) << data.getDiffN();
627 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
628 }
629#endif
630
631 if constexpr (I) {
632
633 auto t_h = getFTensor0FromVec(*hPtr);
634 auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
635
636 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
637
638 const double r = t_coords(0);
639 const double alpha = t_w * vol * cylindrical(r);
640
641 const double set_h = init_h(t_coords(0), t_coords(1), t_coords(2));
642 const double m = get_M(set_h) * alpha;
643
644 int bb = 0;
645 for (; bb != nbRows; ++bb) {
646 locF[bb] += (t_base * alpha) * (t_h - set_h);
647 locF[bb] += (t_diff_base(i) * m) * t_grad_g(i);
648 ++t_base;
649 ++t_diff_base;
650 }
651
652 for (; bb < nbRowBaseFunctions; ++bb) {
653 ++t_base;
654 ++t_diff_base;
655 }
656
657 ++t_h;
658 ++t_grad_g;
659
660 ++t_coords;
661 ++t_w;
662 }
663
664 } else {
665
666 auto t_dot_h = getFTensor0FromVec(*dotHPtr);
667 auto t_h = getFTensor0FromVec(*hPtr);
668 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
669 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
670 auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
671
672 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
673
674 const double r = t_coords(0);
675 const double alpha = t_w * vol * cylindrical(r);
676
677 const double m = get_M(t_h) * alpha;
678
679 int bb = 0;
680 for (; bb != nbRows; ++bb) {
681 locF[bb] += (t_base * alpha) * (t_dot_h);
682 locF[bb] += (t_base * alpha) * (t_grad_h(i) * t_u(i));
683 locF[bb] += (t_diff_base(i) * t_grad_g(i)) * m;
684 ++t_base;
685 ++t_diff_base;
686 }
687
688 for (; bb < nbRowBaseFunctions; ++bb) {
689 ++t_base;
690 ++t_diff_base;
691 }
692
693 ++t_dot_h;
694 ++t_h;
695 ++t_grad_g;
696 ++t_u;
697 ++t_grad_h;
698
699 ++t_coords;
700 ++t_w;
701 }
702 }
703
705 }
706
707private:
708 boost::shared_ptr<MatrixDouble> uPtr;
709 boost::shared_ptr<VectorDouble> dotHPtr;
710 boost::shared_ptr<VectorDouble> hPtr;
711 boost::shared_ptr<MatrixDouble> gradHPtr;
712 boost::shared_ptr<MatrixDouble> gradGPtr;
713};
714
715/**
716 * @brief Lhs for H dU
717 *
718 */
720 OpLhsH_dU(const std::string h_field_name, const std::string u_field_name,
721 boost::shared_ptr<MatrixDouble> grad_h_ptr)
722 : AssemblyDomainEleOp(h_field_name, u_field_name,
723 AssemblyDomainEleOp::OPROWCOL),
724 gradHPtr(grad_h_ptr) {
725 sYmm = false;
726 assembleTranspose = false;
727 }
728 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
729 EntitiesFieldData::EntData &col_data) {
731
732 const double vol = getMeasure();
733 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
734 auto t_coords = getFTensor1CoordsAtGaussPts();
735
736 auto t_row_base = row_data.getFTensor0N();
737 auto t_w = getFTensor0IntegrationWeight();
738
739 for (int gg = 0; gg != nbIntegrationPts; gg++) {
740
741 const double alpha = t_w * vol;
742 auto t_mat = getFTensor1FromPtr<U_FIELD_DIM>(&locMat(0, 0));
743
744 int rr = 0;
745 for (; rr != nbRows; ++rr) {
746 auto t_col_base = col_data.getFTensor0N(gg, 0);
747 for (int cc = 0; cc != nbCols / U_FIELD_DIM; ++cc) {
748 t_mat(i) += (t_row_base * t_col_base * alpha) * t_grad_h(i);
749 ++t_mat;
750 ++t_col_base;
751 }
752 ++t_row_base;
753 }
754
755 for (; rr < nbRowBaseFunctions; ++rr)
756 ++t_row_base;
757
758 ++t_grad_h;
759 ++t_w;
760 }
761
763 }
764
765private:
766 boost::shared_ptr<MatrixDouble> gradHPtr;
767};
768
769/**
770 * @brief Lhs for H dH
771 *
772 */
773template <bool I> struct OpLhsH_dH : public AssemblyDomainEleOp {
774
775 OpLhsH_dH(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
776 boost::shared_ptr<VectorDouble> h_ptr,
777 boost::shared_ptr<MatrixDouble> grad_g_ptr)
779 AssemblyDomainEleOp::OPROWCOL),
780 uPtr(u_ptr), hPtr(h_ptr), gradGPtr(grad_g_ptr) {
781 sYmm = false;
782 }
783
784 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
785 EntitiesFieldData::EntData &col_data) {
787
788 const double vol = getMeasure();
789 auto t_w = getFTensor0IntegrationWeight();
790 auto t_coords = getFTensor1CoordsAtGaussPts();
791 auto t_row_base = row_data.getFTensor0N();
792 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
793
794#ifndef NDEBUG
795 if (row_data.getDiffN().size1() != row_data.getN().size1())
796 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
797 if (row_data.getDiffN().size2() != row_data.getN().size2() * SPACE_DIM) {
798 MOFEM_LOG("SELF", Sev::error)
799 << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
800 MOFEM_LOG("SELF", Sev::error) << row_data.getN();
801 MOFEM_LOG("SELF", Sev::error) << row_data.getDiffN();
802 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
803 }
804
805 if (col_data.getDiffN().size1() != col_data.getN().size1())
806 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
807 if (col_data.getDiffN().size2() != col_data.getN().size2() * SPACE_DIM) {
808 MOFEM_LOG("SELF", Sev::error)
809 << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
810 MOFEM_LOG("SELF", Sev::error) << col_data.getN();
811 MOFEM_LOG("SELF", Sev::error) << col_data.getDiffN();
812 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
813 }
814#endif
815
816 if constexpr (I) {
817
818 auto t_h = getFTensor0FromVec(*hPtr);
819 auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
820
821 for (int gg = 0; gg != nbIntegrationPts; gg++) {
822
823 const double r = t_coords(0);
824 const double alpha = t_w * vol * cylindrical(r);
825
826 int rr = 0;
827 for (; rr != nbRows; ++rr) {
828
829 auto t_col_base = col_data.getFTensor0N(gg, 0);
830 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
831
832 for (int cc = 0; cc != nbCols; ++cc) {
833
834 locMat(rr, cc) += (t_row_base * t_col_base * alpha);
835
836 ++t_col_base;
837 ++t_col_diff_base;
838 }
839
840 ++t_row_base;
841 ++t_row_diff_base;
842 }
843
844 for (; rr < nbRowBaseFunctions; ++rr) {
845 ++t_row_base;
846 ++t_row_diff_base;
847 }
848
849 ++t_h;
850 ++t_grad_g;
851 ++t_w;
852 ++t_coords;
853 }
854
855 } else {
856
857 auto t_h = getFTensor0FromVec(*hPtr);
858 auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
859 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
860
861 auto ts_a = getTSa();
862
863 for (int gg = 0; gg != nbIntegrationPts; gg++) {
864
865 const double r = t_coords(0);
866 const double alpha = t_w * vol * cylindrical(r);
867
868 auto m_dh = get_M_dh(t_h) * alpha;
869 int rr = 0;
870 for (; rr != nbRows; ++rr) {
871
872 auto t_col_base = col_data.getFTensor0N(gg, 0);
873 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
874
875 for (int cc = 0; cc != nbCols; ++cc) {
876
877 locMat(rr, cc) += (t_row_base * t_col_base * alpha) * ts_a;
878 locMat(rr, cc) +=
879 (t_row_base * alpha) * (t_col_diff_base(i) * t_u(i));
880 locMat(rr, cc) +=
881 (t_row_diff_base(i) * t_grad_g(i)) * (t_col_base * m_dh);
882
883 ++t_col_base;
884 ++t_col_diff_base;
885 }
886
887 ++t_row_base;
888 ++t_row_diff_base;
889 }
890
891 for (; rr < nbRowBaseFunctions; ++rr) {
892 ++t_row_base;
893 ++t_row_diff_base;
894 }
895
896 ++t_u;
897 ++t_h;
898 ++t_grad_g;
899 ++t_w;
900 ++t_coords;
901 }
902 }
903
905 }
906
907private:
908 boost::shared_ptr<MatrixDouble> uPtr;
909 boost::shared_ptr<VectorDouble> hPtr;
910 boost::shared_ptr<MatrixDouble> gradGPtr;
911};
912
913/**
914 * @brief Lhs for H dH
915 *
916 */
917template <bool I> struct OpLhsH_dG : public AssemblyDomainEleOp {
918
919 OpLhsH_dG(const std::string field_name_h, const std::string field_name_g,
920 boost::shared_ptr<VectorDouble> h_ptr)
921 : AssemblyDomainEleOp(field_name_h, field_name_g,
922 AssemblyDomainEleOp::OPROWCOL),
923 hPtr(h_ptr) {
924 sYmm = false;
925 assembleTranspose = false;
926 }
927
928 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
929 EntitiesFieldData::EntData &col_data) {
931
932 const double vol = getMeasure();
933 auto t_h = getFTensor0FromVec(*hPtr);
934
935 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
936 auto t_w = getFTensor0IntegrationWeight();
937 auto t_coords = getFTensor1CoordsAtGaussPts();
938
939 for (int gg = 0; gg != nbIntegrationPts; gg++) {
940
941 const double r = t_coords(0);
942 const double alpha = t_w * vol * cylindrical(r);
943
944 double set_h;
945 if constexpr (I)
946 set_h = init_h(t_coords(0), t_coords(1), t_coords(2));
947 else
948 set_h = t_h;
949
950 auto m = get_M(set_h) * alpha;
951
952 int rr = 0;
953 for (; rr != nbRows; ++rr) {
954 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
955
956 for (int cc = 0; cc != nbCols; ++cc) {
957 locMat(rr, cc) += (t_row_diff_base(i) * t_col_diff_base(i)) * m;
958
959 ++t_col_diff_base;
960 }
961
962 ++t_row_diff_base;
963 }
964
965 for (; rr < nbRowBaseFunctions; ++rr) {
966 ++t_row_diff_base;
967 }
968
969 ++t_h;
970 ++t_w;
971 ++t_coords;
972 }
973
975 }
976
977private:
978 boost::shared_ptr<VectorDouble> hPtr;
979};
980
981template <bool I> struct OpRhsG : public AssemblyDomainEleOp {
982
983 OpRhsG(const std::string field_name, boost::shared_ptr<VectorDouble> h_ptr,
984 boost::shared_ptr<MatrixDouble> grad_h_ptr,
985 boost::shared_ptr<VectorDouble> g_ptr)
987 hPtr(h_ptr), gradHPtr(grad_h_ptr), gPtr(g_ptr) {}
988
989 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
991
992 const double vol = getMeasure();
993 auto t_h = getFTensor0FromVec(*hPtr);
994 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
995 auto t_g = getFTensor0FromVec(*gPtr);
996 auto t_coords = getFTensor1CoordsAtGaussPts();
997
998 auto t_base = data.getFTensor0N();
999 auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
1000 auto t_w = getFTensor0IntegrationWeight();
1001
1002 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1003
1004 const double r = t_coords(0);
1005 const double alpha = t_w * vol * cylindrical(r);
1006
1007 double set_h;
1008 if constexpr (I)
1009 set_h = init_h(t_coords(0), t_coords(1), t_coords(2));
1010 else
1011 set_h = t_h;
1012
1013 const double f = get_f(set_h);
1014
1015 int bb = 0;
1016 for (; bb != nbRows; ++bb) {
1017 locF[bb] += (t_base * alpha) * (t_g - f);
1018 locF[bb] -= (t_diff_base(i) * (eta2 * alpha)) * t_grad_h(i);
1019 ++t_base;
1020 ++t_diff_base;
1021 }
1022
1023 for (; bb < nbRowBaseFunctions; ++bb) {
1024 ++t_base;
1025 ++t_diff_base;
1026 }
1027
1028 ++t_h;
1029 ++t_grad_h;
1030 ++t_g;
1031
1032 ++t_coords;
1033 ++t_w;
1034 }
1035
1037 }
1038
1039private:
1040 boost::shared_ptr<VectorDouble> hPtr;
1041 boost::shared_ptr<MatrixDouble> gradHPtr;
1042 boost::shared_ptr<VectorDouble> gPtr;
1043};
1044
1045/**
1046 * @brief Lhs for H dH
1047 *
1048 */
1049template <bool I> struct OpLhsG_dH : public AssemblyDomainEleOp {
1050
1051 OpLhsG_dH(const std::string field_name_g, const std::string field_name_h,
1052 boost::shared_ptr<VectorDouble> h_ptr)
1053 : AssemblyDomainEleOp(field_name_g, field_name_h,
1054 AssemblyDomainEleOp::OPROWCOL),
1055 hPtr(h_ptr) {
1056 sYmm = false;
1057 assembleTranspose = false;
1058 }
1059
1060 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1061 EntitiesFieldData::EntData &col_data) {
1063
1064 const double vol = getMeasure();
1065 auto t_h = getFTensor0FromVec(*hPtr);
1066
1067 auto t_row_base = row_data.getFTensor0N();
1068 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
1069 auto t_w = getFTensor0IntegrationWeight();
1070 auto t_coords = getFTensor1CoordsAtGaussPts();
1071
1072 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1073
1074 const double r = t_coords(0);
1075 const double alpha = t_w * vol * cylindrical(r);
1076
1077 const double f_dh = get_f_dh(t_h) * alpha;
1078 const double beta = eta2 * alpha;
1079
1080 int rr = 0;
1081 for (; rr != nbRows; ++rr) {
1082
1083 auto t_col_base = col_data.getFTensor0N(gg, 0);
1084 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1085
1086 for (int cc = 0; cc != nbCols; ++cc) {
1087
1088 if constexpr (I == false)
1089 locMat(rr, cc) -= (t_row_base * t_col_base) * f_dh;
1090 locMat(rr, cc) -= (t_row_diff_base(i) * beta) * t_col_diff_base(i);
1091
1092 ++t_col_base;
1093 ++t_col_diff_base;
1094 }
1095
1096 ++t_row_base;
1097 ++t_row_diff_base;
1098 }
1099
1100 for (; rr < nbRowBaseFunctions; ++rr) {
1101 ++t_row_base;
1102 ++t_row_diff_base;
1103 }
1104
1105 ++t_h;
1106 ++t_w;
1107 ++t_coords;
1108 }
1109
1111 }
1112
1113private:
1114 boost::shared_ptr<VectorDouble> hPtr;
1115};
1116
1117/**
1118 * @brief Lhs for H dH
1119 *
1120 */
1122
1123 OpLhsG_dG(const std::string field_name)
1125 AssemblyDomainEleOp::OPROWCOL) {
1126 sYmm = true;
1127 }
1128
1129 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1130 EntitiesFieldData::EntData &col_data) {
1132
1133 const double vol = getMeasure();
1134
1135 auto t_row_base = row_data.getFTensor0N();
1136 auto t_w = getFTensor0IntegrationWeight();
1137 auto t_coords = getFTensor1CoordsAtGaussPts();
1138
1139 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1140
1141 const double r = t_coords(0);
1142 const double alpha = t_w * vol * cylindrical(r);
1143
1144 int rr = 0;
1145 for (; rr != nbRows; ++rr) {
1146 auto t_col_base = col_data.getFTensor0N(gg, 0);
1147 const double beta = alpha * t_row_base;
1148 for (int cc = 0; cc != nbCols; ++cc) {
1149 locMat(rr, cc) += (t_col_base * beta);
1150 ++t_col_base;
1151 }
1152
1153 ++t_row_base;
1154 }
1155
1156 for (; rr < nbRowBaseFunctions; ++rr) {
1157 ++t_row_base;
1158 }
1159
1160 ++t_w;
1161 ++t_coords;
1162 }
1163
1165 }
1166
1167private:
1168};
1169
1170
1171
1172} // namespace FreeSurfaceOps
static Index< 'I', 3 > I
constexpr double a
constexpr int SPACE_DIM
Kronecker Delta class symmetric.
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
@ CYLINDRICAL
Definition: definitions.h:117
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
FTensor::Index< 'm', SPACE_DIM > m
auto cylindrical
constexpr int U_FIELD_DIM
constexpr double rho_ave
constexpr CoordinateTypes coord_type
select coordinate system <CARTESIAN, CYLINDRICAL>;
constexpr double mu_ave
auto get_M_dh
constexpr double mu_diff
auto init_h
constexpr auto t_kd
auto d_phase_function_h
auto get_f_dh
auto get_M
auto phase_function
constexpr double rho_diff
auto get_D
constexpr double eta2
const double kappa
auto get_f
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
constexpr double a0
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
double rho
Definition: plastic.cpp:98
constexpr auto field_name
constexpr double mu
boost::shared_ptr< VectorDouble > liftPtr
OpCalculateLift(const std::string field_name, boost::shared_ptr< VectorDouble > p_ptr, boost::shared_ptr< VectorDouble > lift_ptr, boost::shared_ptr< Range > ents_ptr)
boost::shared_ptr< VectorDouble > pPtr
boost::shared_ptr< Range > entsPtr
MoFEMErrorCode doWork(int row_side, EntityType row_type, HookeElement::EntData &data)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpLhsG_dG(const std::string field_name)
boost::shared_ptr< VectorDouble > hPtr
OpLhsG_dH(const std::string field_name_g, const std::string field_name_h, boost::shared_ptr< VectorDouble > h_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpLhsH_dG(const std::string field_name_h, const std::string field_name_g, boost::shared_ptr< VectorDouble > h_ptr)
boost::shared_ptr< VectorDouble > hPtr
boost::shared_ptr< MatrixDouble > uPtr
OpLhsH_dH(const std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< VectorDouble > h_ptr, boost::shared_ptr< MatrixDouble > grad_g_ptr)
boost::shared_ptr< MatrixDouble > gradGPtr
boost::shared_ptr< VectorDouble > hPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< MatrixDouble > gradHPtr
OpLhsH_dU(const std::string h_field_name, const std::string u_field_name, boost::shared_ptr< MatrixDouble > grad_h_ptr)
OpLhsU_dG(const std::string field_name_u, const std::string field_name_h, boost::shared_ptr< MatrixDouble > grad_h_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< MatrixDouble > gradHPtr
boost::shared_ptr< VectorDouble > gPtr
OpLhsU_dH(const std::string field_name_u, const std::string field_name_h, boost::shared_ptr< MatrixDouble > dot_u_ptr, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > grad_u_ptr, boost::shared_ptr< VectorDouble > h_ptr, boost::shared_ptr< VectorDouble > g_ptr)
boost::shared_ptr< VectorDouble > hPtr
boost::shared_ptr< MatrixDouble > uPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< MatrixDouble > gradUPtr
boost::shared_ptr< MatrixDouble > dotUPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
boost::shared_ptr< MatrixDouble > gradUPtr
OpLhsU_dU(const std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > grad_u_ptr, boost::shared_ptr< VectorDouble > h_ptr)
boost::shared_ptr< MatrixDouble > uPtr
boost::shared_ptr< VectorDouble > hPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
OpNormalConstrainLhs(const std::string field_name_row, const std::string field_name_col)
OpNormalConstrainRhs(const std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr)
boost::shared_ptr< MatrixDouble > uPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
boost::shared_ptr< VectorDouble > lambdaPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
OpNormalForceRhs(const std::string field_name, boost::shared_ptr< VectorDouble > lambda_ptr)
OpRhsG(const std::string field_name, boost::shared_ptr< VectorDouble > h_ptr, boost::shared_ptr< MatrixDouble > grad_h_ptr, boost::shared_ptr< VectorDouble > g_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< MatrixDouble > gradHPtr
boost::shared_ptr< VectorDouble > hPtr
boost::shared_ptr< VectorDouble > gPtr
boost::shared_ptr< VectorDouble > hPtr
boost::shared_ptr< MatrixDouble > gradGPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< MatrixDouble > uPtr
boost::shared_ptr< VectorDouble > dotHPtr
boost::shared_ptr< MatrixDouble > gradHPtr
OpRhsH(const std::string field_name, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< VectorDouble > dot_h_ptr, boost::shared_ptr< VectorDouble > h_ptr, boost::shared_ptr< MatrixDouble > grad_h_ptr, boost::shared_ptr< MatrixDouble > grad_g_ptr)
boost::shared_ptr< MatrixDouble > gradUPtr
boost::shared_ptr< VectorDouble > pPtr
OpRhsU(const std::string field_name, boost::shared_ptr< MatrixDouble > dot_u_ptr, boost::shared_ptr< MatrixDouble > u_ptr, boost::shared_ptr< MatrixDouble > grad_u_ptr, boost::shared_ptr< VectorDouble > h_ptr, boost::shared_ptr< MatrixDouble > grad_h_ptr, boost::shared_ptr< VectorDouble > g_ptr, boost::shared_ptr< VectorDouble > p_ptr)
boost::shared_ptr< VectorDouble > gPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data)
boost::shared_ptr< VectorDouble > hPtr
boost::shared_ptr< MatrixDouble > dotUPtr
boost::shared_ptr< MatrixDouble > uPtr
boost::shared_ptr< MatrixDouble > gradHPtr
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
EntityHandle getFEEntityHandle() const
Return finite element entity handle.
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element