v0.14.0
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;
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 OpWettingAngleRhs(const std::string field_name,
140 boost::shared_ptr<MatrixDouble> grad_h_ptr,
141 boost::shared_ptr<Range> ents_ptr = nullptr,
142 double wetting_angle = 0)
144 AssemblyBoundaryEleOp::OPROW),
145 gradHPtr(grad_h_ptr), entsPtr(ents_ptr), wettingAngle(wetting_angle) {}
146
147 MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data) {
149 if (entsPtr) {
150 if (entsPtr->find(AssemblyBoundaryEleOp::getFEEntityHandle()) ==
151 entsPtr->end())
153 }
154 const double area = getMeasure();
155 auto t_w = getFTensor0IntegrationWeight();
156 auto t_row_base = row_data.getFTensor0N();
157 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
158 auto t_coords = getFTensor1CoordsAtGaussPts();
159
160 auto s = wetting_angle_sub_stepping(getFEMethod()->ts_step);
161
162 for (int gg = 0; gg != nbIntegrationPts; gg++) {
163
164 const double r = t_coords(0);
165 const double alpha = t_w * cylindrical(r) * area;
166 const double h_grad_norm = sqrt(t_grad_h(i) * t_grad_h(i) + eps);
167 const double cos_angle = std::cos(M_PI * wettingAngle / 180);
168 const double rhs_wetting = s * eta2 * h_grad_norm * cos_angle;
169
170 // cerr << "pass "
171 // << h_grad_norm <<"\n";
172 int bb = 0;
173 for (; bb != nbRows; ++bb) {
174 locF[bb] += alpha * t_row_base * rhs_wetting;
175 ++t_row_base;
176 }
177
178 for (; bb < nbRowBaseFunctions; ++bb)
179 ++t_row_base;
180
181 ++t_w;
182 ++t_grad_h;
183 ++t_coords;
184 }
185
187 }
188
189private:
190 boost::shared_ptr<MatrixDouble> gradHPtr;
191 boost::shared_ptr<Range> entsPtr;
193};
194
196
197 OpNormalConstrainLhs(const std::string field_name_row,
198 const std::string field_name_col)
199 : AssemblyBoundaryEleOp(field_name_row, field_name_col,
200 AssemblyBoundaryEleOp::OPROWCOL) {
201 assembleTranspose = true;
202 sYmm = false;
203 }
204
205 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
206 EntitiesFieldData::EntData &col_data) {
208
209 auto t_w = getFTensor0IntegrationWeight();
210 auto t_normal = getFTensor1Normal();
211 auto t_row_base = row_data.getFTensor0N();
212 auto t_coords = getFTensor1CoordsAtGaussPts();
213
214 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
215
216 auto t_mat = getFTensor1FromPtr<U_FIELD_DIM>(&locMat(0, 0));
217
218 const double r = t_coords(0);
219 const double alpha = t_w * cylindrical(r);
220
221 int rr = 0;
222 for (; rr != nbRows; ++rr) {
223
224 auto t_col_base = col_data.getFTensor0N(gg, 0);
225 const double a = alpha * t_row_base;
226
227 for (int cc = 0; cc != nbCols / U_FIELD_DIM; ++cc) {
228 t_mat(i) += (a * t_col_base) * t_normal(i);
229 ++t_col_base;
230 ++t_mat;
231 }
232 ++t_row_base;
233 }
234
235 for (; rr < nbRowBaseFunctions; ++rr)
236 ++t_row_base;
237
238 ++t_w;
239 ++t_coords;
240 }
241
243 };
244};
245
247
249 const std::string row_field_name,
250 boost::shared_ptr<MatrixDouble> grad_h_ptr,
251 boost::shared_ptr<std::vector<VectorInt>> col_ind_ptr,
252 boost::shared_ptr<std::vector<MatrixDouble>> col_diff_base_ptr,
253 boost::shared_ptr<Range> ents_ptr = nullptr, double wetting_angle = 0)
254 : BoundaryEleOp(row_field_name, BoundaryEleOp::OPROW),
255 gradHPtr(grad_h_ptr), colIndicesPtr(col_ind_ptr),
256 colDiffBaseFunctionsPtr(col_diff_base_ptr), entsPtr(ents_ptr),
258
259 MoFEMErrorCode doWork(int side, EntityType type,
260 DataForcesAndSourcesCore::EntData &data) {
262 if (entsPtr) {
263 if (entsPtr->find(BoundaryEleOp::getFEEntityHandle()) == entsPtr->end())
265 }
266 const double area = getMeasure();
267
268 const auto row_size = data.getIndices().size();
269 if (row_size == 0)
271
272 auto integrate = [&](auto col_indicies, auto &col_diff_base_functions) {
274
275 const auto col_size = col_indicies.size();
276
277 locMat.resize(row_size, col_size, false);
278 locMat.clear();
279 int nb_gp = getGaussPts().size2();
280 int nb_rows = data.getIndices().size();
281
282 auto t_w = getFTensor0IntegrationWeight();
283 auto t_coords = getFTensor1CoordsAtGaussPts();
284 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
285 auto t_row_base = data.getFTensor0N();
286 int nb_row_base_functions = data.getN().size2();
287
288 auto s = wetting_angle_sub_stepping(getFEMethod()->ts_step);
289
290 for (int gg = 0; gg != nb_gp; ++gg) {
291
292 const double r = t_coords(0);
293 const double alpha = t_w * area * cylindrical(r);
294 const double h_grad_norm = sqrt(t_grad_h(i) * t_grad_h(i) + eps);
295 const double one_over_h_grad_norm = 1. / h_grad_norm;
296 const double beta = s * alpha * eta2 * one_over_h_grad_norm *
297 std::cos(M_PI * wettingAngle / 180);
298
299 int rr = 0;
300 for (; rr != nb_rows; ++rr) {
301 const double delta = beta * t_row_base;
302
303 auto ptr = &col_diff_base_functions(gg, 0);
304 auto t_col_diff_base = getFTensor1FromPtr<SPACE_DIM>(ptr);
305
306 for (int cc = 0; cc != col_size; ++cc) {
307 locMat(rr, cc) += t_col_diff_base(i) * (delta * t_grad_h(i));
308 ++t_col_diff_base;
309 }
310 ++t_row_base;
311 }
312
313 for (; rr < nb_row_base_functions; ++rr) {
314 ++t_row_base;
315 }
316
317 ++t_grad_h;
318 ++t_w;
319 ++t_coords;
320 }
321
323 };
324
325 for (auto c = 0; c != colIndicesPtr->size(); ++c) {
326
327 auto &col_ind = (*colIndicesPtr)[c];
328 if (col_ind.size()) {
329 auto &diff_base = (*colDiffBaseFunctionsPtr)[c];
330
331 CHKERR integrate(col_ind, diff_base);
332
333 CHKERR MatSetValues(getKSPB(), data.getIndices().size(),
334 &*data.getIndices().begin(), col_ind.size(),
335 &*col_ind.begin(), &locMat(0, 0), ADD_VALUES);
336 }
337 }
338
340 }
341
342private:
343 MatrixDouble locMat;
344
345 boost::shared_ptr<MatrixDouble> gradHPtr;
346 boost::shared_ptr<Range> entsPtr;
348 boost::shared_ptr<std::vector<VectorInt>> colIndicesPtr;
349 boost::shared_ptr<std::vector<MatrixDouble>> colDiffBaseFunctionsPtr;
350};
351
352/**
353 * @brief Rhs for U
354 *
355 */
357
358 OpRhsU(const std::string field_name,
359 boost::shared_ptr<MatrixDouble> dot_u_ptr,
360 boost::shared_ptr<MatrixDouble> u_ptr,
361 boost::shared_ptr<MatrixDouble> grad_u_ptr,
362 boost::shared_ptr<VectorDouble> h_ptr,
363 boost::shared_ptr<MatrixDouble> grad_h_ptr,
364 boost::shared_ptr<VectorDouble> g_ptr,
365 boost::shared_ptr<VectorDouble> p_ptr)
367 dotUPtr(dot_u_ptr), uPtr(u_ptr), gradUPtr(grad_u_ptr), hPtr(h_ptr),
368 gradHPtr(grad_h_ptr), gPtr(g_ptr), pPtr(p_ptr) {}
369
370 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
372
373 const double vol = getMeasure();
374 auto t_dot_u = getFTensor1FromMat<U_FIELD_DIM>(*dotUPtr);
375 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
376 auto t_p = getFTensor0FromVec(*pPtr);
377 auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*gradUPtr);
378 auto t_h = getFTensor0FromVec(*hPtr);
379 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
380 auto t_g = getFTensor0FromVec(*gPtr);
381 auto t_coords = getFTensor1CoordsAtGaussPts();
382
383 auto t_base = data.getFTensor0N();
384 auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
385
386 auto t_w = getFTensor0IntegrationWeight();
387
396
397 t_buoyancy(i) = 0;
398 t_gravity(i) = 0;
399
400 for (int gg = 0; gg != nbIntegrationPts; gg++) {
401
402 const double r = t_coords(0);
403 const double alpha = t_w * vol * cylindrical(r);
404
405 const double rho = phase_function(t_h, rho_diff, rho_ave);
406 const double mu = phase_function(t_h, mu_diff, mu_ave);
407
408 auto t_D = get_D(2 * mu);
409
410 t_inertia_force(i) = (rho * alpha) * (t_dot_u(i));
411 // t_buoyancy(SPACE_DIM - 1) = -(alpha * rho * a0) * t_h;
412 t_gravity(SPACE_DIM - 1) = (alpha * rho * a0);
413 t_phase_force(i) = -alpha * kappa * t_g * t_grad_h(i);
414 t_convection(i) = (rho * alpha) * (t_u(j) * t_grad_u(i, j));
415
416 t_stress(i, j) =
417 alpha * (t_D(i, j, k, l) * t_grad_u(k, l) + t_kd(i, j) * t_p);
418
419 auto t_nf = getFTensor1FromArray<U_FIELD_DIM, U_FIELD_DIM>(locF);
420
421 t_forces(i) = t_inertia_force(i) + t_buoyancy(i) + t_gravity(i) +
422 t_convection(i) + t_phase_force(i);
423
424 int bb = 0;
425 for (; bb != nbRows / U_FIELD_DIM; ++bb) {
426
427 t_nf(i) += t_base * t_forces(i);
428 t_nf(i) += t_diff_base(j) * t_stress(i, j);
429
430 if (coord_type == CYLINDRICAL) {
431 t_nf(0) += (t_base * (alpha / t_coords(0))) * (2 * mu * t_u(0) + t_p);
432 }
433
434 ++t_base;
435 ++t_diff_base;
436 ++t_nf;
437 }
438
439 for (; bb < nbRowBaseFunctions; ++bb) {
440 ++t_diff_base;
441 ++t_base;
442 }
443
444 ++t_dot_u;
445 ++t_u;
446 ++t_grad_u;
447 ++t_h;
448 ++t_grad_h;
449 ++t_g;
450 ++t_p;
451
452 ++t_w;
453 ++t_coords;
454 }
455
457 }
458
459private:
460 boost::shared_ptr<MatrixDouble> dotUPtr;
461 boost::shared_ptr<MatrixDouble> uPtr;
462 boost::shared_ptr<MatrixDouble> gradUPtr;
463 boost::shared_ptr<VectorDouble> hPtr;
464 boost::shared_ptr<MatrixDouble> gradHPtr;
465 boost::shared_ptr<VectorDouble> gPtr;
466 boost::shared_ptr<VectorDouble> pPtr;
467};
468
469/**
470 * @brief Lhs for U dU
471 *
472 */
474
475 OpLhsU_dU(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
476 boost::shared_ptr<MatrixDouble> grad_u_ptr,
477 boost::shared_ptr<VectorDouble> h_ptr)
479 AssemblyDomainEleOp::OPROWCOL),
480 uPtr(u_ptr), gradUPtr(grad_u_ptr), hPtr(h_ptr) {
481 sYmm = false;
482 assembleTranspose = false;
483 }
484
485 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
486 EntitiesFieldData::EntData &col_data) {
488
489 const double vol = getMeasure();
490 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
491 auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*gradUPtr);
492 auto t_h = getFTensor0FromVec(*hPtr);
493 auto t_coords = getFTensor1CoordsAtGaussPts();
494
495 auto t_row_base = row_data.getFTensor0N();
496 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
497
498 auto t_w = getFTensor0IntegrationWeight();
499
500 auto get_mat = [&](const int rr) {
501 return getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(locMat, rr);
502 };
503
504 auto ts_a = getTSa();
506
507 for (int gg = 0; gg != nbIntegrationPts; gg++) {
508
509 const double r = t_coords(0);
510 const double alpha = t_w * vol * cylindrical(r);
511 const double rho = phase_function(t_h, rho_diff, rho_ave);
512 const double mu = phase_function(t_h, mu_diff, mu_ave);
513
514 const double beta0 = alpha * rho;
515 const double beta1 = beta0 * ts_a;
516 auto t_D = get_D(2 * mu);
517
518 int rr = 0;
519 for (; rr != nbRows / U_FIELD_DIM; ++rr) {
520
521 auto t_mat = get_mat(rr * U_FIELD_DIM);
522 auto t_col_base = col_data.getFTensor0N(gg, 0);
523 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
524
526 // I mix up the indices here so that it behaves like a
527 // Dg. That way I don't have to have a separate wrapper
528 // class Christof_Expr, which simplifies things.
529 t_d_stress(l, j, k) = t_D(i, j, k, l) * (alpha * t_row_diff_base(i));
530
531 for (int cc = 0; cc != nbCols / U_FIELD_DIM; ++cc) {
532
533 const double bb = t_row_base * t_col_base;
534
535 t_mat(i, j) += (beta1 * bb) * t_kd(i, j);
536 t_mat(i, j) += (beta0 * bb) * t_grad_u(i, j);
537 t_mat(i, j) +=
538 (beta0 * t_row_base) * t_kd(i, j) * (t_col_diff_base(k) * t_u(k));
539 t_mat(i, j) += t_d_stress(i, j, k) * t_col_diff_base(k);
540
541 if (coord_type == CYLINDRICAL) {
542 t_mat(0, 0) += (bb * (alpha / t_coords(0))) * (2 * mu);
543 }
544
545 ++t_mat;
546 ++t_col_base;
547 ++t_col_diff_base;
548 }
549
550 ++t_row_base;
551 ++t_row_diff_base;
552 }
553
554 for (; rr < nbRowBaseFunctions; ++rr) {
555 ++t_row_diff_base;
556 ++t_row_base;
557 }
558
559 ++t_u;
560 ++t_grad_u;
561 ++t_h;
562
563 ++t_coords;
564 ++t_w;
565 }
566
568 }
569
570private:
571 boost::shared_ptr<MatrixDouble> uPtr;
572 boost::shared_ptr<MatrixDouble> gradUPtr;
573 boost::shared_ptr<VectorDouble> hPtr;
574};
575
576struct OpLoopSideGetDataForSideEle : ForcesAndSourcesCore::UserDataOperator {
577
578 using UDO = ForcesAndSourcesCore::UserDataOperator;
579
581 const std::string field_name,
582 boost::shared_ptr<std::vector<VectorInt>> col_indices_ptr,
583 boost::shared_ptr<std::vector<MatrixDouble>> col_diff_basefunctions_ptr)
584 : UDO(field_name, UDO::OPCOL), colIndicesPtr(col_indices_ptr),
585 colDiffBaseFunctionsPtr(col_diff_basefunctions_ptr) {}
586
587 MoFEMErrorCode doWork(int side, EntityType type,
588 DataForcesAndSourcesCore::EntData &data) {
590
591 if (type == MBVERTEX) {
592 colIndicesPtr->clear();
594 }
595
596 colIndicesPtr->push_back(data.getIndices());
597 colDiffBaseFunctionsPtr->push_back(data.getDiffN());
598
600 }
601
602protected:
603 boost::shared_ptr<std::vector<VectorInt>> colIndicesPtr;
604 boost::shared_ptr<std::vector<MatrixDouble>> colDiffBaseFunctionsPtr;
605};
606
607/**
608 * @brief Lhs for U dH
609 *
610 */
612
613 OpLhsU_dH(const std::string field_name_u, const std::string field_name_h,
614 boost::shared_ptr<MatrixDouble> dot_u_ptr,
615 boost::shared_ptr<MatrixDouble> u_ptr,
616 boost::shared_ptr<MatrixDouble> grad_u_ptr,
617 boost::shared_ptr<VectorDouble> h_ptr,
618 boost::shared_ptr<VectorDouble> g_ptr)
619 : AssemblyDomainEleOp(field_name_u, field_name_h,
620 AssemblyDomainEleOp::OPROWCOL),
621 dotUPtr(dot_u_ptr), uPtr(u_ptr), gradUPtr(grad_u_ptr), hPtr(h_ptr),
622 gPtr(g_ptr) {
623 sYmm = false;
624 assembleTranspose = false;
625 }
626
627 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
628 EntitiesFieldData::EntData &col_data) {
630
631 const double vol = getMeasure();
632 auto t_dot_u = getFTensor1FromMat<U_FIELD_DIM>(*dotUPtr);
633 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
634 auto t_grad_u = getFTensor2FromMat<U_FIELD_DIM, SPACE_DIM>(*gradUPtr);
635 auto t_h = getFTensor0FromVec(*hPtr);
636 auto t_g = getFTensor0FromVec(*gPtr);
637 auto t_coords = getFTensor1CoordsAtGaussPts();
638
639 auto t_row_base = row_data.getFTensor0N();
640 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
641
642 auto t_w = getFTensor0IntegrationWeight();
643
646 FTensor::Tensor1<double, U_FIELD_DIM> t_inertia_force_dh;
651
652 t_buoyancy_dh(i) = 0;
653 t_gravity_dh(i) = 0;
654
655 for (int gg = 0; gg != nbIntegrationPts; gg++) {
656
657 const double r = t_coords(0);
658 const double alpha = t_w * vol * cylindrical(r);
659
660 const double rho = phase_function(t_h, rho_diff, rho_ave);
661 const double rho_dh = d_phase_function_h(t_h, rho_diff);
662 const double mu_dh = d_phase_function_h(t_h, mu_diff);
663
664 auto t_D_dh = get_D(2 * mu_dh);
665
666 t_inertia_force_dh(i) = (alpha * rho_dh) * t_dot_u(i);
667 // t_buoyancy_dh(SPACE_DIM - 1) = -(alpha * a0) * (rho + rho_dh * t_h);
668 t_gravity_dh(SPACE_DIM - 1) = (alpha * rho_dh * a0);
669 t_convection_dh(i) = (rho_dh * alpha) * (t_u(j) * t_grad_u(i, j));
670 const double t_phase_force_g_dh = -alpha * kappa * t_g;
671 t_forces_dh(i) = t_inertia_force_dh(i) + t_buoyancy_dh(i) +
672 t_gravity_dh(i) + t_convection_dh(i);
673
674 t_stress_dh(i, j) = alpha * (t_D_dh(i, j, k, l) * t_grad_u(k, l));
675
676 int rr = 0;
677 for (; rr != nbRows / U_FIELD_DIM; ++rr) {
678
679 auto t_mat =
680 getFTensor1FromMat<U_FIELD_DIM, 1>(locMat, rr * U_FIELD_DIM);
681 auto t_col_base = col_data.getFTensor0N(gg, 0);
682 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
683
684 for (int cc = 0; cc != nbCols; ++cc) {
685
686 const double bb = t_row_base * t_col_base;
687 t_mat(i) += t_forces_dh(i) * bb;
688 t_mat(i) += (t_phase_force_g_dh * t_row_base) * t_col_diff_base(i);
689 t_mat(i) += (t_row_diff_base(j) * t_col_base) * t_stress_dh(i, j);
690
691 if (coord_type == CYLINDRICAL) {
692 t_mat(0) += (bb * (alpha / t_coords(0))) * (2 * mu_dh * t_u(0));
693 }
694
695 ++t_mat;
696 ++t_col_base;
697 ++t_col_diff_base;
698 }
699
700 ++t_row_base;
701 ++t_row_diff_base;
702 }
703
704 for (; rr < nbRowBaseFunctions; ++rr) {
705 ++t_row_diff_base;
706 ++t_row_base;
707 }
708
709 ++t_dot_u;
710 ++t_u;
711 ++t_grad_u;
712 ++t_h;
713 ++t_g;
714 ++t_coords;
715 ++t_w;
716 }
717
719 }
720
721private:
722 boost::shared_ptr<MatrixDouble> dotUPtr;
723 boost::shared_ptr<MatrixDouble> uPtr;
724 boost::shared_ptr<MatrixDouble> gradUPtr;
725 boost::shared_ptr<VectorDouble> hPtr;
726 boost::shared_ptr<VectorDouble> gPtr;
727};
728
729/**
730 * @brief Lhs for G dH
731 *
732 */
734
735 OpLhsU_dG(const std::string field_name_u, const std::string field_name_h,
736 boost::shared_ptr<MatrixDouble> grad_h_ptr)
737 : AssemblyDomainEleOp(field_name_u, field_name_h,
738 AssemblyDomainEleOp::OPROWCOL),
739 gradHPtr(grad_h_ptr) {
740 sYmm = false;
741 assembleTranspose = false;
742 }
743
744 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
745 EntitiesFieldData::EntData &col_data) {
747
748 const double vol = getMeasure();
749 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
750 auto t_coords = getFTensor1CoordsAtGaussPts();
751
752 auto t_row_base = row_data.getFTensor0N();
753 auto t_w = getFTensor0IntegrationWeight();
754
755 for (int gg = 0; gg != nbIntegrationPts; gg++) {
756
757 const double r = t_coords(0);
758 const double alpha = t_w * vol * cylindrical(r);
759
761 t_phase_force_dg(i) = -alpha * kappa * t_grad_h(i);
762
763 int rr = 0;
764 for (; rr != nbRows / U_FIELD_DIM; ++rr) {
765 auto t_mat =
766 getFTensor1FromMat<U_FIELD_DIM, 1>(locMat, rr * U_FIELD_DIM);
767 auto t_col_base = col_data.getFTensor0N(gg, 0);
768
769 for (int cc = 0; cc != nbCols; ++cc) {
770 const double bb = t_row_base * t_col_base;
771 t_mat(i) += t_phase_force_dg(i) * bb;
772
773 ++t_mat;
774 ++t_col_base;
775 }
776
777 ++t_row_base;
778 }
779
780 for (; rr < nbRowBaseFunctions; ++rr) {
781 ++t_row_base;
782 }
783
784 ++t_grad_h;
785 ++t_coords;
786 ++t_w;
787 }
788
790 }
791
792private:
793 boost::shared_ptr<MatrixDouble> gradHPtr;
794};
795
796template <bool I> struct OpRhsH : public AssemblyDomainEleOp {
797
798 OpRhsH(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
799 boost::shared_ptr<VectorDouble> dot_h_ptr,
800 boost::shared_ptr<VectorDouble> h_ptr,
801 boost::shared_ptr<MatrixDouble> grad_h_ptr,
802 boost::shared_ptr<MatrixDouble> grad_g_ptr)
804 uPtr(u_ptr), dotHPtr(dot_h_ptr), hPtr(h_ptr), gradHPtr(grad_h_ptr),
805 gradGPtr(grad_g_ptr) {}
806
807 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
809
810 const double vol = getMeasure();
811 auto t_w = getFTensor0IntegrationWeight();
812 auto t_coords = getFTensor1CoordsAtGaussPts();
813 auto t_base = data.getFTensor0N();
814 auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
815
816#ifndef NDEBUG
817 if (data.getDiffN().size1() != data.getN().size1())
818 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
819 if (data.getDiffN().size2() != data.getN().size2() * SPACE_DIM) {
820 MOFEM_LOG("SELF", Sev::error)
821 << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
822 MOFEM_LOG("SELF", Sev::error) << data.getN();
823 MOFEM_LOG("SELF", Sev::error) << data.getDiffN();
824 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
825 }
826#endif
827
828 if constexpr (I) {
829
830 auto t_h = getFTensor0FromVec(*hPtr);
831 auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
832
833 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
834
835 const double r = t_coords(0);
836 const double alpha = t_w * vol * cylindrical(r);
837
838 const double set_h = init_h(t_coords(0), t_coords(1), t_coords(2));
839 const double m = get_M(set_h) * alpha;
840
841 int bb = 0;
842 for (; bb != nbRows; ++bb) {
843 locF[bb] += (t_base * alpha) * (t_h - set_h);
844 locF[bb] += (t_diff_base(i) * m) * t_grad_g(i);
845 ++t_base;
846 ++t_diff_base;
847 }
848
849 for (; bb < nbRowBaseFunctions; ++bb) {
850 ++t_base;
851 ++t_diff_base;
852 }
853
854 ++t_h;
855 ++t_grad_g;
856
857 ++t_coords;
858 ++t_w;
859 }
860
861 } else {
862
863 auto t_dot_h = getFTensor0FromVec(*dotHPtr);
864 auto t_h = getFTensor0FromVec(*hPtr);
865 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
866 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
867 auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
868
869 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
870
871 const double r = t_coords(0);
872 const double alpha = t_w * vol * cylindrical(r);
873
874 const double m = get_M(t_h) * alpha;
875
876 int bb = 0;
877 for (; bb != nbRows; ++bb) {
878 locF[bb] += (t_base * alpha) * (t_dot_h);
879 locF[bb] += (t_base * alpha) * (t_grad_h(i) * t_u(i));
880 locF[bb] += (t_diff_base(i) * t_grad_g(i)) * m;
881 ++t_base;
882 ++t_diff_base;
883 }
884
885 for (; bb < nbRowBaseFunctions; ++bb) {
886 ++t_base;
887 ++t_diff_base;
888 }
889
890 ++t_dot_h;
891 ++t_h;
892 ++t_grad_g;
893 ++t_u;
894 ++t_grad_h;
895
896 ++t_coords;
897 ++t_w;
898 }
899 }
900
902 }
903
904private:
905 boost::shared_ptr<MatrixDouble> uPtr;
906 boost::shared_ptr<VectorDouble> dotHPtr;
907 boost::shared_ptr<VectorDouble> hPtr;
908 boost::shared_ptr<MatrixDouble> gradHPtr;
909 boost::shared_ptr<MatrixDouble> gradGPtr;
910};
911
912/**
913 * @brief Lhs for H dU
914 *
915 */
917 OpLhsH_dU(const std::string h_field_name, const std::string u_field_name,
918 boost::shared_ptr<MatrixDouble> grad_h_ptr)
919 : AssemblyDomainEleOp(h_field_name, u_field_name,
920 AssemblyDomainEleOp::OPROWCOL),
921 gradHPtr(grad_h_ptr) {
922 sYmm = false;
923 assembleTranspose = false;
924 }
925 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
926 EntitiesFieldData::EntData &col_data) {
928
929 const double vol = getMeasure();
930 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
931 auto t_coords = getFTensor1CoordsAtGaussPts();
932
933 auto t_row_base = row_data.getFTensor0N();
934 auto t_w = getFTensor0IntegrationWeight();
935
936 for (int gg = 0; gg != nbIntegrationPts; gg++) {
937
938 const auto r = t_coords(0);
939 const auto alpha = t_w * vol * cylindrical(r);
940 auto t_mat = getFTensor1FromPtr<U_FIELD_DIM>(&locMat(0, 0));
942
943 int rr = 0;
944 for (; rr != nbRows; ++rr) {
945 t_row(i) = (alpha * t_row_base) * t_grad_h(i);
946 auto t_col_base = col_data.getFTensor0N(gg, 0);
947 for (int cc = 0; cc != nbCols / U_FIELD_DIM; ++cc) {
948 t_mat(i) += t_row(i) * t_col_base;
949 ++t_mat;
950 ++t_col_base;
951 }
952 ++t_row_base;
953 }
954
955 for (; rr < nbRowBaseFunctions; ++rr)
956 ++t_row_base;
957
958 ++t_grad_h;
959 ++t_w;
960 ++t_coords;
961 }
962
964 }
965
966private:
967 boost::shared_ptr<MatrixDouble> gradHPtr;
968};
969
970/**
971 * @brief Lhs for H dH
972 *
973 */
974template <bool I> struct OpLhsH_dH : public AssemblyDomainEleOp {
975
976 OpLhsH_dH(const std::string field_name, boost::shared_ptr<MatrixDouble> u_ptr,
977 boost::shared_ptr<VectorDouble> h_ptr,
978 boost::shared_ptr<MatrixDouble> grad_g_ptr)
980 AssemblyDomainEleOp::OPROWCOL),
981 uPtr(u_ptr), hPtr(h_ptr), gradGPtr(grad_g_ptr) {
982 sYmm = false;
983 }
984
985 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
986 EntitiesFieldData::EntData &col_data) {
988
989 const double vol = getMeasure();
990 auto t_w = getFTensor0IntegrationWeight();
991 auto t_coords = getFTensor1CoordsAtGaussPts();
992 auto t_row_base = row_data.getFTensor0N();
993 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
994
995#ifndef NDEBUG
996 if (row_data.getDiffN().size1() != row_data.getN().size1())
997 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
998 if (row_data.getDiffN().size2() != row_data.getN().size2() * SPACE_DIM) {
999 MOFEM_LOG("SELF", Sev::error)
1000 << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
1001 MOFEM_LOG("SELF", Sev::error) << row_data.getN();
1002 MOFEM_LOG("SELF", Sev::error) << row_data.getDiffN();
1003 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
1004 }
1005
1006 if (col_data.getDiffN().size1() != col_data.getN().size1())
1007 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 1");
1008 if (col_data.getDiffN().size2() != col_data.getN().size2() * SPACE_DIM) {
1009 MOFEM_LOG("SELF", Sev::error)
1010 << "Side " << rowSide << " " << CN::EntityTypeName(rowType);
1011 MOFEM_LOG("SELF", Sev::error) << col_data.getN();
1012 MOFEM_LOG("SELF", Sev::error) << col_data.getDiffN();
1013 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "wrong size 2");
1014 }
1015#endif
1016
1017 if constexpr (I) {
1018
1019 auto t_h = getFTensor0FromVec(*hPtr);
1020 auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
1021
1022 for (int gg = 0; gg != nbIntegrationPts; gg++) {
1023
1024 const double r = t_coords(0);
1025 const double alpha = t_w * vol * cylindrical(r);
1026
1027 int rr = 0;
1028 for (; rr != nbRows; ++rr) {
1029
1030 auto t_col_base = col_data.getFTensor0N(gg, 0);
1031 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1032
1033 for (int cc = 0; cc != nbCols; ++cc) {
1034
1035 locMat(rr, cc) += (t_row_base * t_col_base * alpha);
1036
1037 ++t_col_base;
1038 ++t_col_diff_base;
1039 }
1040
1041 ++t_row_base;
1042 ++t_row_diff_base;
1043 }
1044
1045 for (; rr < nbRowBaseFunctions; ++rr) {
1046 ++t_row_base;
1047 ++t_row_diff_base;
1048 }
1049
1050 ++t_h;
1051 ++t_grad_g;
1052 ++t_w;
1053 ++t_coords;
1054 }
1055
1056 } else {
1057
1058 auto t_h = getFTensor0FromVec(*hPtr);
1059 auto t_grad_g = getFTensor1FromMat<SPACE_DIM>(*gradGPtr);
1060 auto t_u = getFTensor1FromMat<U_FIELD_DIM>(*uPtr);
1061
1062 auto ts_a = getTSa();
1063
1064 for (int gg = 0; gg != nbIntegrationPts; gg++) {
1065
1066 const double r = t_coords(0);
1067 const double alpha = t_w * vol * cylindrical(r);
1068
1069 auto m_dh = get_M_dh(t_h) * alpha;
1070 int rr = 0;
1071 for (; rr != nbRows; ++rr) {
1072
1073 auto t_col_base = col_data.getFTensor0N(gg, 0);
1074 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1075
1076 for (int cc = 0; cc != nbCols; ++cc) {
1077
1078 locMat(rr, cc) += (t_row_base * t_col_base * alpha) * ts_a;
1079 locMat(rr, cc) +=
1080 (t_row_base * alpha) * (t_col_diff_base(i) * t_u(i));
1081 locMat(rr, cc) +=
1082 (t_row_diff_base(i) * t_grad_g(i)) * (t_col_base * m_dh);
1083
1084 ++t_col_base;
1085 ++t_col_diff_base;
1086 }
1087
1088 ++t_row_base;
1089 ++t_row_diff_base;
1090 }
1091
1092 for (; rr < nbRowBaseFunctions; ++rr) {
1093 ++t_row_base;
1094 ++t_row_diff_base;
1095 }
1096
1097 ++t_u;
1098 ++t_h;
1099 ++t_grad_g;
1100 ++t_w;
1101 ++t_coords;
1102 }
1103 }
1104
1106 }
1107
1108private:
1109 boost::shared_ptr<MatrixDouble> uPtr;
1110 boost::shared_ptr<VectorDouble> hPtr;
1111 boost::shared_ptr<MatrixDouble> gradGPtr;
1112};
1113
1114/**
1115 * @brief Lhs for H dH
1116 *
1117 */
1118template <bool I> struct OpLhsH_dG : public AssemblyDomainEleOp {
1119
1120 OpLhsH_dG(const std::string field_name_h, const std::string field_name_g,
1121 boost::shared_ptr<VectorDouble> h_ptr)
1122 : AssemblyDomainEleOp(field_name_h, field_name_g,
1123 AssemblyDomainEleOp::OPROWCOL),
1124 hPtr(h_ptr) {
1125 sYmm = false;
1126 assembleTranspose = false;
1127 }
1128
1129 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1130 EntitiesFieldData::EntData &col_data) {
1132
1133 const double vol = getMeasure();
1134 auto t_h = getFTensor0FromVec(*hPtr);
1135
1136 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
1137 auto t_w = getFTensor0IntegrationWeight();
1138 auto t_coords = getFTensor1CoordsAtGaussPts();
1139
1140 for (int gg = 0; gg != nbIntegrationPts; gg++) {
1141
1142 const double r = t_coords(0);
1143 const double alpha = t_w * vol * cylindrical(r);
1144
1145 double set_h;
1146 if constexpr (I)
1147 set_h = init_h(t_coords(0), t_coords(1), t_coords(2));
1148 else
1149 set_h = t_h;
1150
1151 auto m = get_M(set_h) * alpha;
1152
1153 int rr = 0;
1154 for (; rr != nbRows; ++rr) {
1155 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1156
1157 for (int cc = 0; cc != nbCols; ++cc) {
1158 locMat(rr, cc) += (t_row_diff_base(i) * t_col_diff_base(i)) * m;
1159
1160 ++t_col_diff_base;
1161 }
1162
1163 ++t_row_diff_base;
1164 }
1165
1166 for (; rr < nbRowBaseFunctions; ++rr) {
1167 ++t_row_diff_base;
1168 }
1169
1170 ++t_h;
1171 ++t_w;
1172 ++t_coords;
1173 }
1174
1176 }
1177
1178private:
1179 boost::shared_ptr<VectorDouble> hPtr;
1180};
1181
1182template <bool I> struct OpRhsG : public AssemblyDomainEleOp {
1183
1184 OpRhsG(const std::string field_name, boost::shared_ptr<VectorDouble> h_ptr,
1185 boost::shared_ptr<MatrixDouble> grad_h_ptr,
1186 boost::shared_ptr<VectorDouble> g_ptr)
1188 hPtr(h_ptr), gradHPtr(grad_h_ptr), gPtr(g_ptr) {}
1189
1190 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &data) {
1192
1193 const double vol = getMeasure();
1194 auto t_h = getFTensor0FromVec(*hPtr);
1195 auto t_grad_h = getFTensor1FromMat<SPACE_DIM>(*gradHPtr);
1196 auto t_g = getFTensor0FromVec(*gPtr);
1197 auto t_coords = getFTensor1CoordsAtGaussPts();
1198
1199 auto t_base = data.getFTensor0N();
1200 auto t_diff_base = data.getFTensor1DiffN<SPACE_DIM>();
1201 auto t_w = getFTensor0IntegrationWeight();
1202
1203 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1204
1205 const double r = t_coords(0);
1206 const double alpha = t_w * vol * cylindrical(r);
1207
1208 double set_h;
1209 if constexpr (I)
1210 set_h = init_h(t_coords(0), t_coords(1), t_coords(2));
1211 else
1212 set_h = t_h;
1213
1214 const double f = get_f(set_h);
1215
1216 int bb = 0;
1217 for (; bb != nbRows; ++bb) {
1218 locF[bb] += (t_base * alpha) * (t_g - f);
1219 locF[bb] -= (t_diff_base(i) * (eta2 * alpha)) * t_grad_h(i);
1220 ++t_base;
1221 ++t_diff_base;
1222 }
1223
1224 for (; bb < nbRowBaseFunctions; ++bb) {
1225 ++t_base;
1226 ++t_diff_base;
1227 }
1228
1229 ++t_h;
1230 ++t_grad_h;
1231 ++t_g;
1232
1233 ++t_coords;
1234 ++t_w;
1235 }
1236
1238 }
1239
1240private:
1241 boost::shared_ptr<VectorDouble> hPtr;
1242 boost::shared_ptr<MatrixDouble> gradHPtr;
1243 boost::shared_ptr<VectorDouble> gPtr;
1244};
1245
1246/**
1247 * @brief Lhs for H dH
1248 *
1249 */
1250template <bool I> struct OpLhsG_dH : public AssemblyDomainEleOp {
1251
1252 OpLhsG_dH(const std::string field_name_g, const std::string field_name_h,
1253 boost::shared_ptr<VectorDouble> h_ptr)
1254 : AssemblyDomainEleOp(field_name_g, field_name_h,
1255 AssemblyDomainEleOp::OPROWCOL),
1256 hPtr(h_ptr) {
1257 sYmm = false;
1258 assembleTranspose = false;
1259 }
1260
1261 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1262 EntitiesFieldData::EntData &col_data) {
1264
1265 const double vol = getMeasure();
1266 auto t_h = getFTensor0FromVec(*hPtr);
1267
1268 auto t_row_base = row_data.getFTensor0N();
1269 auto t_row_diff_base = row_data.getFTensor1DiffN<SPACE_DIM>();
1270 auto t_w = getFTensor0IntegrationWeight();
1271 auto t_coords = getFTensor1CoordsAtGaussPts();
1272
1273 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1274
1275 const double r = t_coords(0);
1276 const double alpha = t_w * vol * cylindrical(r);
1277
1278 const double f_dh = get_f_dh(t_h) * alpha;
1279 const double beta = eta2 * alpha;
1280
1281 int rr = 0;
1282 for (; rr != nbRows; ++rr) {
1283
1284 auto t_col_base = col_data.getFTensor0N(gg, 0);
1285 auto t_col_diff_base = col_data.getFTensor1DiffN<SPACE_DIM>(gg, 0);
1286
1287 for (int cc = 0; cc != nbCols; ++cc) {
1288
1289 if constexpr (I == false)
1290 locMat(rr, cc) -= (t_row_base * t_col_base) * f_dh;
1291 locMat(rr, cc) -= (t_row_diff_base(i) * beta) * t_col_diff_base(i);
1292
1293 ++t_col_base;
1294 ++t_col_diff_base;
1295 }
1296
1297 ++t_row_base;
1298 ++t_row_diff_base;
1299 }
1300
1301 for (; rr < nbRowBaseFunctions; ++rr) {
1302 ++t_row_base;
1303 ++t_row_diff_base;
1304 }
1305
1306 ++t_h;
1307 ++t_w;
1308 ++t_coords;
1309 }
1310
1312 }
1313
1314private:
1315 boost::shared_ptr<VectorDouble> hPtr;
1316};
1317
1318/**
1319 * @brief Lhs for H dH
1320 *
1321 */
1323
1324 OpLhsG_dG(const std::string field_name)
1326 AssemblyDomainEleOp::OPROWCOL) {
1327 sYmm = true;
1328 }
1329
1330 MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data,
1331 EntitiesFieldData::EntData &col_data) {
1333
1334 const double vol = getMeasure();
1335
1336 auto t_row_base = row_data.getFTensor0N();
1337 auto t_w = getFTensor0IntegrationWeight();
1338 auto t_coords = getFTensor1CoordsAtGaussPts();
1339
1340 for (int gg = 0; gg != nbIntegrationPts; ++gg) {
1341
1342 const double r = t_coords(0);
1343 const double alpha = t_w * vol * cylindrical(r);
1344
1345 int rr = 0;
1346 for (; rr != nbRows; ++rr) {
1347 auto t_col_base = col_data.getFTensor0N(gg, 0);
1348 const double beta = alpha * t_row_base;
1349 for (int cc = 0; cc != nbCols; ++cc) {
1350 locMat(rr, cc) += (t_col_base * beta);
1351 ++t_col_base;
1352 }
1353
1354 ++t_row_base;
1355 }
1356
1357 for (; rr < nbRowBaseFunctions; ++rr) {
1358 ++t_row_base;
1359 }
1360
1361 ++t_w;
1362 ++t_coords;
1363 }
1364
1366 }
1367
1368private:
1369};
1370
1371} // namespace FreeSurfaceOps
constexpr double a
static const double eps
constexpr int SPACE_DIM
Kronecker Delta class symmetric.
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#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
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'm', SPACE_DIM > m
auto cylindrical
double kappa
double mu_diff
constexpr int U_FIELD_DIM
auto get_M_dh
double rho_diff
auto init_h
Initialisation function.
constexpr auto t_kd
auto d_phase_function_h
double rho_ave
double eta2
auto get_f_dh
auto wetting_angle
auto get_M
double mu_ave
auto phase_function
auto get_D
auto wetting_angle_sub_stepping
auto get_f
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
constexpr double a0
FTensor::Index< 'i', SPACE_DIM > i
constexpr CoordinateTypes coord_type
const double c
speed of light (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
constexpr IntegrationType I
double rho
Definition: plastic.cpp:191
constexpr auto field_name
static constexpr double delta
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
boost::shared_ptr< std::vector< VectorInt > > colIndicesPtr
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
OpLoopSideGetDataForSideEle(const std::string field_name, boost::shared_ptr< std::vector< VectorInt > > col_indices_ptr, boost::shared_ptr< std::vector< MatrixDouble > > col_diff_basefunctions_ptr)
ForcesAndSourcesCore::UserDataOperator UDO
boost::shared_ptr< std::vector< MatrixDouble > > colDiffBaseFunctionsPtr
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
boost::shared_ptr< MatrixDouble > gradHPtr
boost::shared_ptr< std::vector< VectorInt > > colIndicesPtr
boost::shared_ptr< std::vector< MatrixDouble > > colDiffBaseFunctionsPtr
boost::shared_ptr< Range > entsPtr
OpWettingAngleLhs(const std::string row_field_name, boost::shared_ptr< MatrixDouble > grad_h_ptr, boost::shared_ptr< std::vector< VectorInt > > col_ind_ptr, boost::shared_ptr< std::vector< MatrixDouble > > col_diff_base_ptr, boost::shared_ptr< Range > ents_ptr=nullptr, double wetting_angle=0)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
boost::shared_ptr< MatrixDouble > gradHPtr
OpWettingAngleRhs(const std::string field_name, boost::shared_ptr< MatrixDouble > grad_h_ptr, boost::shared_ptr< Range > ents_ptr=nullptr, double wetting_angle=0)
MoFEMErrorCode iNtegrate(DataForcesAndSourcesCore::EntData &row_data)
boost::shared_ptr< Range > entsPtr
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.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
const FEMethod * getFEMethod() const
Return raw pointer to Finite Element Method object.
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element