v0.14.0
Loading...
Searching...
No Matches
NavierStokesElement.cpp
Go to the documentation of this file.
1/**
2 * \file NavierStokesElement.cpp
3 * \example NavierStokesElement.cpp
4 *
5 * \brief Implementation of operators for fluid flow
6 *
7 * FIXME: Code does not handle higher order geometry
8 *
9 * Implementation of operators for simulation of the fluid flow governed by
10 * Stokes and Navier-Stokes equations, and computation of the drag
11 * force on a fluid-solid inteface
12 */
13
14
15
16#include <MoFEM.hpp>
17using namespace MoFEM;
18using namespace boost::numeric;
21
23 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_rhs_ptr,
24 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs_ptr,
25 const std::string velocity_field, const std::string pressure_field,
26 boost::shared_ptr<CommonData> common_data, const EntityType type) {
28
29 for (auto &sit : common_data->setOfBlocksData) {
30
31 if (type == MBPRISM) {
32 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
33 fe_lhs_ptr->getOpPtrVector().push_back(
35 fe_lhs_ptr->getOpPtrVector().push_back(
36 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
37 fe_lhs_ptr->getOpPtrVector().push_back(
38 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
39 fe_rhs_ptr->getOpPtrVector().push_back(
41 fe_rhs_ptr->getOpPtrVector().push_back(
42 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
43 fe_rhs_ptr->getOpPtrVector().push_back(
44 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
45 }
46
47 fe_lhs_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
48 velocity_field, common_data->velPtr));
49 fe_lhs_ptr->getOpPtrVector().push_back(
50 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
51 common_data->gradVelPtr));
52
53 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsDiagLin(
54 velocity_field, velocity_field, common_data, sit.second));
55 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsDiagNonLin(
56 velocity_field, velocity_field, common_data, sit.second));
57 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsOffDiag(
58 velocity_field, pressure_field, common_data, sit.second));
59
60 fe_rhs_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
61 velocity_field, common_data->velPtr));
62 fe_rhs_ptr->getOpPtrVector().push_back(
63 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
64 common_data->gradVelPtr));
65 fe_rhs_ptr->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
66 pressure_field, common_data->pressPtr));
67
68 fe_rhs_ptr->getOpPtrVector().push_back(
69 new OpAssembleRhsVelocityLin(velocity_field, common_data, sit.second));
70 fe_rhs_ptr->getOpPtrVector().push_back(new OpAssembleRhsVelocityNonLin(
71 velocity_field, common_data, sit.second));
72 fe_rhs_ptr->getOpPtrVector().push_back(
73 new OpAssembleRhsPressure(pressure_field, common_data, sit.second));
74 }
75
77};
78
80 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_rhs_ptr,
81 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs_ptr,
82 const std::string velocity_field, const std::string pressure_field,
83 boost::shared_ptr<CommonData> common_data, const EntityType type) {
85
86 for (auto &sit : common_data->setOfBlocksData) {
87
88 if (type == MBPRISM) {
89 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
90 fe_lhs_ptr->getOpPtrVector().push_back(
92 fe_lhs_ptr->getOpPtrVector().push_back(
93 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
94 fe_lhs_ptr->getOpPtrVector().push_back(
95 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
96 fe_rhs_ptr->getOpPtrVector().push_back(
98 fe_rhs_ptr->getOpPtrVector().push_back(
99 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
100 fe_rhs_ptr->getOpPtrVector().push_back(
101 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
102 }
103
104 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsDiagLin(
105 velocity_field, velocity_field, common_data, sit.second));
106 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsOffDiag(
107 velocity_field, pressure_field, common_data, sit.second));
108
109 fe_rhs_ptr->getOpPtrVector().push_back(
110 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
111 common_data->gradVelPtr));
112
113 fe_rhs_ptr->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
114 pressure_field, common_data->pressPtr));
115 fe_rhs_ptr->getOpPtrVector().push_back(
116 new OpAssembleRhsVelocityLin(velocity_field, common_data, sit.second));
117 fe_rhs_ptr->getOpPtrVector().push_back(
118 new OpAssembleRhsPressure(pressure_field, common_data, sit.second));
119 }
120
122};
123
125 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_flux_ptr,
126 const std::string velocity_field, boost::shared_ptr<CommonData> common_data,
127 const EntityType type) {
129
130 for (auto &sit : common_data->setOfBlocksData) {
131
132 if (type == MBPRISM) {
133 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
134 fe_flux_ptr->getOpPtrVector().push_back(
136 fe_flux_ptr->getOpPtrVector().push_back(
137 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
138 fe_flux_ptr->getOpPtrVector().push_back(
139 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
140 }
141 fe_flux_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
142 velocity_field, common_data->velPtr));
143 fe_flux_ptr->getOpPtrVector().push_back(
144 new OpCalcVolumeFlux(velocity_field, common_data, sit.second));
145 }
146
148};
149
151 boost::shared_ptr<FaceElementForcesAndSourcesCore> dragFe,
152 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideDragFe,
153 std::string side_fe_name, const std::string velocity_field,
154 const std::string pressure_field,
155 boost::shared_ptr<CommonData> common_data) {
157
158 auto det_ptr = boost::make_shared<VectorDouble>();
159 auto jac_ptr = boost::make_shared<MatrixDouble>();
160 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
161
162 for (auto &sit : common_data->setOfFacesData) {
163 sideDragFe->getOpPtrVector().push_back(
164 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
165 common_data->gradVelPtr));
166 dragFe->getOpPtrVector().push_back(new OpCalculateHOJac<2>(jac_ptr));
167 dragFe->getOpPtrVector().push_back(
168 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
169 dragFe->getOpPtrVector().push_back(
170 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
171
172 dragFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
173 pressure_field, common_data->pressPtr));
174
175 dragFe->getOpPtrVector().push_back(
177 velocity_field, sideDragFe, side_fe_name, common_data, sit.second));
178 dragFe->getOpPtrVector().push_back(new NavierStokesElement::OpCalcDragForce(
179 velocity_field, common_data, sit.second));
180 }
182};
183
185 boost::shared_ptr<PostProcFace> postProcDragPtr,
186 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideDragFe,
187 std::string side_fe_name, const std::string velocity_field,
188 const std::string pressure_field,
189 boost::shared_ptr<CommonData> common_data) {
191
192 auto det_ptr = boost::make_shared<VectorDouble>();
193 auto jac_ptr = boost::make_shared<MatrixDouble>();
194 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
195
196 for (auto &sit : common_data->setOfFacesData) {
197 sideDragFe->getOpPtrVector().push_back(
198 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
199 common_data->gradVelPtr));
200 postProcDragPtr->getOpPtrVector().push_back(
201 new OpCalculateHOJac<2>(jac_ptr));
202 postProcDragPtr->getOpPtrVector().push_back(
203 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
204 postProcDragPtr->getOpPtrVector().push_back(
205 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
206
207 postProcDragPtr->getOpPtrVector().push_back(
208 new OpCalculateVectorFieldValues<3>(velocity_field,
209 common_data->velPtr));
210 postProcDragPtr->getOpPtrVector().push_back(
211 new OpCalculateScalarFieldValues(pressure_field,
212 common_data->pressPtr));
213
214 postProcDragPtr->getOpPtrVector().push_back(
216 velocity_field, sideDragFe, side_fe_name, common_data, sit.second));
217 postProcDragPtr->getOpPtrVector().push_back(
219 velocity_field, postProcDragPtr->getPostProcMesh(),
220 postProcDragPtr->getMapGaussPts(), common_data, sit.second));
221 }
223};
224
226 int row_side, int col_side, EntityType row_type, EntityType col_type,
227 EntData &row_data, EntData &col_data) {
228
230
231 row_nb_dofs = row_data.getIndices().size();
232 if (!row_nb_dofs)
234 col_nb_dofs = col_data.getIndices().size();
235 if (!col_nb_dofs)
237
238 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
239 blockData.eNts.end()) {
241 }
242
243 row_nb_gauss_pts = row_data.getN().size1();
244 if (!row_nb_gauss_pts) {
246 }
247
248 isOnDiagonal = (row_type == col_type) && (row_side == col_side);
249
250 // Set size can clear local tangent matrix
251 locMat.resize(row_nb_dofs, col_nb_dofs, false);
252 locMat.clear();
253
254 CHKERR iNtegrate(row_data, col_data);
255
256 CHKERR aSsemble(row_data, col_data);
257
259}
260
262 EntData &col_data) {
264
265 const int *row_ind = &*row_data.getIndices().begin();
266 const int *col_ind = &*col_data.getIndices().begin();
267
268 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
269 : getFEMethod()->snes_B;
270
271 CHKERR MatSetValues(B, row_nb_dofs, row_ind, col_nb_dofs, col_ind,
272 &*locMat.data().begin(), ADD_VALUES);
273
274 if (!diagonalBlock || (sYmm && !isOnDiagonal)) {
275 locMat = trans(locMat);
276 CHKERR MatSetValues(B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
277 &*locMat.data().begin(), ADD_VALUES);
278 }
279
281}
282
285 EntData &col_data) {
287
288 const int row_nb_base_functions = row_data.getN().size2();
289 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
290
292 FTensor::Index<'i', 3> i;
293 FTensor::Index<'j', 3> j;
294
295 // INTEGRATION
296 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
297
298 // Get volume and integration weight
299 double w = getVolume() * getGaussPts()(3, gg);
300
301 int row_bb = 0;
302 for (; row_bb != row_nb_dofs / 3; row_bb++) {
303
304 t1(i) = w * row_diff_base_functions(i);
305
306 auto base_functions = col_data.getFTensor0N(gg, 0);
307 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
308
309 FTensor::Tensor1<double *, 3> k(&locMat(3 * row_bb + 0, col_bb),
310 &locMat(3 * row_bb + 1, col_bb),
311 &locMat(3 * row_bb + 2, col_bb));
312
313 k(i) -= t1(i) * base_functions;
314 ++base_functions;
315 }
316 ++row_diff_base_functions;
317 }
318 for (; row_bb != row_nb_base_functions; row_bb++) {
319 ++row_diff_base_functions;
320 }
321 }
322
324}
325
328 EntData &col_data) {
330
331 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
333 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
334 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
335 &m(r + 2, c + 2));
336 };
337
338 const int row_nb_base_functions = row_data.getN().size2();
339
340 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
341
342 FTensor::Index<'i', 3> i;
343 FTensor::Index<'j', 3> j;
344
345 // integrate local matrix for entity block
346 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
347
348 // Get volume and integration weight
349 double w = getVolume() * getGaussPts()(3, gg);
350 double const alpha = w * blockData.viscousCoef;
351
352 int row_bb = 0;
353 for (; row_bb != row_nb_dofs / 3; row_bb++) {
354
355 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
356
357 const int final_bb = isOnDiagonal ? row_bb + 1 : col_nb_dofs / 3;
358
359 int col_bb = 0;
360
361 for (; col_bb != final_bb; col_bb++) {
362
363 auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
364
365 for (int i : {0, 1, 2}) {
366 for (int j : {0, 1, 2}) {
367 t_assemble(i, i) +=
368 alpha * row_diff_base_functions(j) * col_diff_base_functions(j);
369 t_assemble(i, j) +=
370 alpha * row_diff_base_functions(j) * col_diff_base_functions(i);
371 }
372 }
373
374 // Next base function for column
375 ++col_diff_base_functions;
376 }
377
378 // Next base function for row
379 ++row_diff_base_functions;
380 }
381
382 for (; row_bb != row_nb_base_functions; row_bb++) {
383 ++row_diff_base_functions;
384 }
385 }
386
387 if (isOnDiagonal) {
388 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
389 int col_bb = 0;
390 for (; col_bb != row_bb + 1; col_bb++) {
391 auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
392 auto t_off_side = get_tensor2(locMat, 3 * col_bb, 3 * row_bb);
393 t_off_side(i, j) = t_assemble(j, i);
394 }
395 }
396 }
397
399}
400
403 EntData &col_data) {
405
406 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
408 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
409 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
410 &m(r + 2, c + 2));
411 };
412
413 const int row_nb_base_functions = row_data.getN().size2();
414
415 auto row_base_functions = row_data.getFTensor0N();
416
417 auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
418 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
419
420 FTensor::Index<'i', 3> i;
421 FTensor::Index<'j', 3> j;
422
423 // integrate local matrix for entity block
424 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
425
426 // Get volume and integration weight
427 double w = getVolume() * getGaussPts()(3, gg);
428 double const beta = w * blockData.inertiaCoef;
429
430 int row_bb = 0;
431 for (; row_bb != row_nb_dofs / 3; row_bb++) {
432
433 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
434 auto col_base_functions = col_data.getFTensor0N(gg, 0);
435
436 const int final_bb = col_nb_dofs / 3;
437 int col_bb = 0;
438
439 for (; col_bb != final_bb; col_bb++) {
440
441 auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
442
443 for (int i : {0, 1, 2}) {
444 for (int j : {0, 1, 2}) {
445 t_assemble(i, j) +=
446 beta * col_base_functions * t_u_grad(i, j) * row_base_functions;
447 t_assemble(i, i) +=
448 beta * t_u(j) * col_diff_base_functions(j) * row_base_functions;
449 }
450 }
451
452 // Next base function for column
453 ++col_diff_base_functions;
454 ++col_base_functions;
455 }
456
457 // Next base function for row
458 ++row_base_functions;
459 }
460
461 for (; row_bb != row_nb_base_functions; row_bb++) {
462 ++row_base_functions;
463 }
464
465 ++t_u;
466 ++t_u_grad;
467 }
468
470}
471
473 EntityType row_type,
474 EntData &row_data) {
476 // get number of dofs on row
477 nbRows = row_data.getIndices().size();
478 if (!nbRows)
480 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
481 blockData.eNts.end()) {
483 }
484 // get number of integration points
485 nbIntegrationPts = getGaussPts().size2();
486
487 // integrate local vector
488 CHKERR iNtegrate(row_data);
489 // assemble local vector
490 CHKERR aSsemble(row_data);
492}
493
496 // get global indices of local vector
497 const int *indices = &*data.getIndices().data().begin();
498 // get values from local vector
499 const double *vals = &*locVec.data().begin();
500 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
501 : getFEMethod()->snes_f;
502 // assemble vector
503 CHKERR VecSetValues(f, nbRows, indices, vals, ADD_VALUES);
504
506}
507
511
512 auto get_tensor1 = [](VectorDouble &v, const int r) {
514 &v(r + 0), &v(r + 1), &v(r + 2));
515 };
516
517 // set size of local vector
518 locVec.resize(nbRows, false);
519 // clear local entity vector
520 locVec.clear();
521
522 int nb_base_functions = data.getN().size2();
523
524 // get base function gradient on rows
525 auto t_v_grad = data.getFTensor1DiffN<3>();
526
527 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
528 auto t_p = getFTensor0FromVec(*commonData->pressPtr);
529
530 FTensor::Index<'i', 3> i;
531 FTensor::Index<'j', 3> j;
532
533 // loop over all integration points
534 for (int gg = 0; gg != nbIntegrationPts; gg++) {
535
536 double w = getVolume() * getGaussPts()(3, gg);
537
538 // evaluate constant term
539 const double alpha = w * blockData.viscousCoef;
540
541 auto t_a = get_tensor1(locVec, 0);
542 int rr = 0;
543
544 // loop over base functions
545 for (; rr != nbRows / 3; rr++) {
546 // add to local vector source term
547
548 t_a(i) += alpha * t_u_grad(i, j) * t_v_grad(j);
549 t_a(j) += alpha * t_u_grad(i, j) * t_v_grad(i);
550
551 t_a(i) -= w * t_p * t_v_grad(i);
552
553 ++t_a; // move to next element of local vector
554 ++t_v_grad; // move to next gradient of base function
555 }
556
557 for (; rr != nb_base_functions; rr++) {
558 ++t_v_grad;
559 }
560
561 ++t_u_grad;
562 ++t_p;
563 }
564
566}
567
571
572 auto get_tensor1 = [](VectorDouble &v, const int r) {
574 &v(r + 0), &v(r + 1), &v(r + 2));
575 };
576
577 // set size of local vector
578 locVec.resize(nbRows, false);
579 // clear local entity vector
580 locVec.clear();
581
582 // get base functions on entity
583 auto t_v = data.getFTensor0N();
584
585 int nb_base_functions = data.getN().size2();
586
587 auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
588 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
589
590 FTensor::Index<'i', 3> i;
591 FTensor::Index<'j', 3> j;
592
593 // loop over all integration points
594 for (int gg = 0; gg != nbIntegrationPts; gg++) {
595
596 double w = getVolume() * getGaussPts()(3, gg);
597 // evaluate constant term
598 const double beta = w * blockData.inertiaCoef;
599
600 auto t_a = get_tensor1(locVec, 0);
601 int rr = 0;
602
603 // loop over base functions
604 for (; rr != nbRows / 3; rr++) {
605 // add to local vector source term
606
607 t_a(i) += beta * t_v * t_u_grad(i, j) * t_u(j);
608
609 ++t_a; // move to next element of local vector
610 ++t_v; // move to next base function
611 }
612
613 for (; rr != nb_base_functions; rr++) {
614 ++t_v;
615 }
616
617 ++t_u;
618 ++t_u_grad;
619 }
621}
622
626
627 // commonData->getBlockData(blockData);
628
629 // set size to local vector
630 locVec.resize(nbRows, false);
631 // clear local vector
632 locVec.clear();
633 // get base function
634 auto t_q = data.getFTensor0N();
635
636 int nb_base_functions = data.getN().size2();
637 // get solution at integration point
638
639 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
640
641 // FTensor::Index<'i', 3> i;
642
643 // make loop over integration points
644 for (int gg = 0; gg != nbIntegrationPts; gg++) {
645 // evaluate function on boundary and scale it by area and integration
646 // weight
647
648 double w = getVolume() * getGaussPts()(3, gg);
649
650 // get element of vector
652 int rr = 0;
653 for (; rr != nbRows; rr++) {
654
655 for (int ii : {0, 1, 2}) {
656 t_a -= w * t_q * t_u_grad(ii, ii);
657 }
658
659 ++t_a;
660 ++t_q;
661 }
662
663 for (; rr != nb_base_functions; rr++) {
664 ++t_q;
665 }
666
667 ++t_u_grad;
668 }
670}
671
673 EntityType type,
674 EntData &data) {
676
677 if (type != MBVERTEX)
678 PetscFunctionReturn(0);
679
680 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
681 blockData.eNts.end()) {
683 }
684
685 const int nb_gauss_pts = getGaussPts().size2();
686
687 auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
688 auto t_w = getFTensor0IntegrationWeight(); ///< Integration weight
689
690 FTensor::Index<'i', 3> i;
691
693 t_flux(i) = 0.0; // Zero entries
694
695 // loop over all integration points
696 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
697
698 double vol = getVolume();
699 t_flux(i) += t_w * vol * t_u(i);
700
701 ++t_w;
702 ++t_u;
703 }
704
705 // Set array of indices
706 constexpr std::array<int, 3> indices = {0, 1, 2};
707
708 // Assemble volumetric flux
709 CHKERR VecSetValues(commonData->volumeFluxVec, 3, indices.data(), &t_flux(0),
710 ADD_VALUES);
711
713}
714
716 EntityType type,
717 EntData &data) {
719 if (type != MBVERTEX)
720 PetscFunctionReturn(0);
721
722 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
723 blockData.eNts.end()) {
725 }
726
727 const int nb_gauss_pts = getGaussPts().size2();
728
729 auto pressure_drag_at_gauss_pts =
730 getFTensor1FromMat<3>(*commonData->pressureDragTract);
731 auto shear_drag_at_gauss_pts =
732 getFTensor1FromMat<3>(*commonData->shearDragTract);
733 auto total_drag_at_gauss_pts =
734 getFTensor1FromMat<3>(*commonData->totalDragTract);
735
736 FTensor::Tensor1<double, 3> t_pressure_drag_force;
737 FTensor::Tensor1<double, 3> t_shear_drag_force;
738 FTensor::Tensor1<double, 3> t_total_drag_force;
739
740 FTensor::Index<'i', 3> i;
741
742 t_pressure_drag_force(i) = 0.0; // Zero entries
743 t_shear_drag_force(i) = 0.0; // Zero entries
744 t_total_drag_force(i) = 0.0; // Zero entries
745
746 auto t_w = getFTensor0IntegrationWeight();
747 auto t_normal = getFTensor1NormalsAtGaussPts();
748
749 for (int gg = 0; gg != nb_gauss_pts; gg++) {
750
751 double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
752 double w = t_w * nrm2 * 0.5;
753
754 t_pressure_drag_force(i) += w * pressure_drag_at_gauss_pts(i);
755 t_shear_drag_force(i) += w * shear_drag_at_gauss_pts(i);
756 t_total_drag_force(i) += w * total_drag_at_gauss_pts(i);
757
758 ++t_w;
759 ++t_normal;
760 ++pressure_drag_at_gauss_pts;
761 ++shear_drag_at_gauss_pts;
762 ++total_drag_at_gauss_pts;
763 }
764
765 // Set array of indices
766 constexpr std::array<int, 3> indices = {0, 1, 2};
767
768 CHKERR VecSetValues(commonData->pressureDragForceVec, 3, indices.data(),
769 &t_pressure_drag_force(0), ADD_VALUES);
770 CHKERR VecSetValues(commonData->shearDragForceVec, 3, indices.data(),
771 &t_shear_drag_force(0), ADD_VALUES);
772 CHKERR VecSetValues(commonData->totalDragForceVec, 3, indices.data(),
773 &t_total_drag_force(0), ADD_VALUES);
774
776}
777
779 EntityType type,
780 EntData &data) {
782
783 if (type != MBVERTEX)
784 PetscFunctionReturn(0);
785
786 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
787 blockData.eNts.end()) {
789 }
790
791 CHKERR loopSideVolumes(sideFeName, *sideFe);
792
793 const int nb_gauss_pts = getGaussPts().size2();
794
795 auto t_p = getFTensor0FromVec(*commonData->pressPtr);
796 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
797
798 auto t_normal = getFTensor1NormalsAtGaussPts();
799
800 FTensor::Tensor1<double, 3> t_unit_normal;
801
802 FTensor::Index<'i', 3> i;
803 FTensor::Index<'j', 3> j;
804
805 commonData->pressureDragTract->resize(3, nb_gauss_pts, false);
806 commonData->pressureDragTract->clear();
807
808 commonData->shearDragTract->resize(3, nb_gauss_pts, false);
809 commonData->shearDragTract->clear();
810
811 commonData->totalDragTract->resize(3, nb_gauss_pts, false);
812 commonData->totalDragTract->clear();
813
814 auto pressure_drag_at_gauss_pts =
815 getFTensor1FromMat<3>(*commonData->pressureDragTract);
816 auto shear_drag_at_gauss_pts =
817 getFTensor1FromMat<3>(*commonData->shearDragTract);
818 auto total_drag_at_gauss_pts =
819 getFTensor1FromMat<3>(*commonData->totalDragTract);
820
821 for (int gg = 0; gg != nb_gauss_pts; gg++) {
822
823 double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
824 t_unit_normal(i) = t_normal(i) / nrm2;
825
826 double mu = blockData.fluidViscosity;
827
828 pressure_drag_at_gauss_pts(i) = t_p * t_unit_normal(i);
829 shear_drag_at_gauss_pts(i) =
830 -mu * (t_u_grad(i, j) + t_u_grad(j, i)) * t_unit_normal(j);
831 total_drag_at_gauss_pts(i) =
832 pressure_drag_at_gauss_pts(i) + shear_drag_at_gauss_pts(i);
833
834 ++pressure_drag_at_gauss_pts;
835 ++shear_drag_at_gauss_pts;
836 ++total_drag_at_gauss_pts;
837 ++t_p;
838 ++t_u_grad;
839 ++t_normal;
840 }
842}
843
845 EntityType type,
846 EntData &data) {
848 if (type != MBVERTEX)
849 PetscFunctionReturn(0);
850
851 double def_VAL[9];
852 bzero(def_VAL, 9 * sizeof(double));
853
854 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
855 blockData.eNts.end()) {
857 }
858
859 Tag th_velocity;
860 Tag th_pressure;
861 Tag th_velocity_grad;
862 Tag th_shear_drag;
863 Tag th_pressure_drag;
864 Tag th_total_drag;
865
866 CHKERR postProcMesh.tag_get_handle("VELOCITY", 3, MB_TYPE_DOUBLE, th_velocity,
867 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
868 CHKERR postProcMesh.tag_get_handle("PRESSURE", 1, MB_TYPE_DOUBLE, th_pressure,
869 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
870 CHKERR postProcMesh.tag_get_handle("VELOCITY_GRAD", 9, MB_TYPE_DOUBLE,
871 th_velocity_grad,
872 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
873
874 CHKERR postProcMesh.tag_get_handle("PRESSURE_DRAG", 3, MB_TYPE_DOUBLE,
875 th_pressure_drag,
876 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
877 CHKERR postProcMesh.tag_get_handle("SHEAR_DRAG", 3, MB_TYPE_DOUBLE,
878 th_shear_drag,
879 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
880 CHKERR postProcMesh.tag_get_handle("TOTAL_DRAG", 3, MB_TYPE_DOUBLE,
881 th_total_drag,
882 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
883
884 auto t_p = getFTensor0FromVec(*commonData->pressPtr);
885
886 const int nb_gauss_pts = commonData->pressureDragTract->size2();
887
888 MatrixDouble velGradMat;
889 velGradMat.resize(3, 3);
890 VectorDouble velVec;
891 velVec.resize(3);
892 VectorDouble pressDragVec;
893 pressDragVec.resize(3);
894 VectorDouble viscDragVec;
895 viscDragVec.resize(3);
896 VectorDouble totDragVec;
897 totDragVec.resize(3);
898
899 for (int gg = 0; gg != nb_gauss_pts; gg++) {
900
901 for (int i : {0, 1, 2}) {
902 for (int j : {0, 1, 2}) {
903 velGradMat(i, j) = (*commonData->gradVelPtr)(i * 3 + j, gg);
904 }
905 velVec(i) = (*commonData->velPtr)(i, gg);
906 pressDragVec(i) = (*commonData->pressureDragTract)(i, gg);
907 viscDragVec(i) = (*commonData->shearDragTract)(i, gg);
908 totDragVec(i) = (*commonData->totalDragTract)(i, gg);
909 }
910 CHKERR postProcMesh.tag_set_data(th_velocity, &mapGaussPts[gg], 1,
911 &velVec(0));
912 CHKERR postProcMesh.tag_set_data(th_pressure, &mapGaussPts[gg], 1, &t_p);
913 CHKERR postProcMesh.tag_set_data(th_velocity_grad, &mapGaussPts[gg], 1,
914 &velGradMat(0, 0));
915
916 CHKERR postProcMesh.tag_set_data(th_pressure_drag, &mapGaussPts[gg], 1,
917 &pressDragVec(0));
918
919 CHKERR postProcMesh.tag_set_data(th_shear_drag, &mapGaussPts[gg], 1,
920 &viscDragVec(0));
921
922 CHKERR postProcMesh.tag_set_data(th_total_drag, &mapGaussPts[gg], 1,
923 &totDragVec(0));
924
925 ++t_p;
926 }
927
929}
930
932 EntityType type,
933 EntData &data) {
935 if (type != MBVERTEX)
936 PetscFunctionReturn(0);
937 double def_VAL[9];
938 bzero(def_VAL, 9 * sizeof(double));
939
940 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
941 blockData.eNts.end()) {
943 }
944 // commonData->getBlockData(blockData);
945
946 Tag th_vorticity;
947 Tag th_q;
948 Tag th_l2;
949 CHKERR postProcMesh.tag_get_handle("VORTICITY", 3, MB_TYPE_DOUBLE,
950 th_vorticity, MB_TAG_CREAT | MB_TAG_SPARSE,
951 def_VAL);
952 CHKERR postProcMesh.tag_get_handle("Q", 1, MB_TYPE_DOUBLE, th_q,
953 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
954 CHKERR postProcMesh.tag_get_handle("L2", 1, MB_TYPE_DOUBLE, th_l2,
955 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
956
957 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
958 // auto p = getFTensor0FromVec(*commonData->pressPtr);
959
960 const int nb_gauss_pts = commonData->gradVelPtr->size2();
961 // const int nb_gauss_pts2 = commonData->pressPtr->size();
962
963 // const double lambda = commonData->lAmbda;
964 // const double mu = commonData->mU;
965
966 // FTensor::Index<'i', 3> i;
967 // FTensor::Index<'j', 3> j;
968 // FTensor::Index<'j', 3> k;
970 // FTensor::Tensor2<double,3,3> t_s;
971 double q;
972 double l2;
973 // FTensor::Tensor2<double, 3, 3> stress;
974 MatrixDouble S;
975 MatrixDouble Omega;
977
978 S.resize(3, 3);
979 Omega.resize(3, 3);
980 M.resize(3, 3);
981
982 for (int gg = 0; gg != nb_gauss_pts; gg++) {
983
984 vorticity(0) = t_u_grad(2, 1) - t_u_grad(1, 2);
985 vorticity(1) = t_u_grad(0, 2) - t_u_grad(2, 0);
986 vorticity(2) = t_u_grad(1, 0) - t_u_grad(0, 1);
987
988 q = 0;
989 for (int i = 0; i != 3; i++) {
990 for (int j = 0; j != 3; j++) {
991 q += -0.5 * t_u_grad(i, j) * t_u_grad(j, i);
992 }
993 }
994 for (int i = 0; i != 3; i++) {
995 for (int j = 0; j != 3; j++) {
996 S(i, j) = 0.5 * (t_u_grad(i, j) + t_u_grad(j, i));
997 Omega(i, j) = 0.5 * (t_u_grad(i, j) - t_u_grad(j, i));
998 M(i, j) = 0.0;
999 }
1000 }
1001
1002 for (int i = 0; i != 3; i++) {
1003 for (int j = 0; j != 3; j++) {
1004 for (int k = 0; k != 3; k++) {
1005 M(i, j) += S(i, k) * S(k, j) + Omega(i, k) * Omega(k, j);
1006 }
1007 }
1008 }
1009
1010 MatrixDouble eigen_vectors = M;
1011 VectorDouble eigen_values(3);
1012
1013 // LAPACK - eigenvalues and vectors. Applied twice for initial creates
1014 // memory space
1015 int n = 3, lda = 3, info, lwork = -1;
1016 double wkopt;
1017 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
1018 &(eigen_values.data()[0]), &wkopt, lwork);
1019 if (info != 0)
1020 SETERRQ1(PETSC_COMM_SELF, 1,
1021 "is something wrong with lapack_dsyev info = %d", info);
1022 lwork = (int)wkopt;
1023 double work[lwork];
1024 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
1025 &(eigen_values.data()[0]), work, lwork);
1026 if (info != 0)
1027 SETERRQ1(PETSC_COMM_SELF, 1,
1028 "is something wrong with lapack_dsyev info = %d", info);
1029
1030 map<double, int> eigen_sort;
1031 eigen_sort[eigen_values[0]] = 0;
1032 eigen_sort[eigen_values[1]] = 1;
1033 eigen_sort[eigen_values[2]] = 2;
1034
1035 // prin_stress_vect.clear();
1036 VectorDouble prin_vals_vect(3);
1037 prin_vals_vect.clear();
1038
1039 int ii = 0;
1040 for (map<double, int>::reverse_iterator mit = eigen_sort.rbegin();
1041 mit != eigen_sort.rend(); mit++) {
1042 prin_vals_vect[ii] = eigen_values[mit->second];
1043 // for (int dd = 0; dd != 3; dd++) {
1044 // prin_stress_vect(ii, dd) = eigen_vectors.data()[3 * mit->second +
1045 // dd];
1046 // }
1047 ii++;
1048 }
1049
1050 l2 = prin_vals_vect[1];
1051 // cout << prin_vals_vect << endl;
1052 // cout << "-0.5 sum: " << -0.5 * (prin_vals_vect[0] + prin_vals_vect[1]
1053 // + prin_vals_vect[2]) << endl; cout << "q: " << q << endl;
1054
1055 // t_s(i,j) = 0.5*(t)
1056
1057 // vorticity(0) = t_u_grad(1, 2) - t_u_grad(2, 1);
1058 // vorticity(1) = t_u_grad(2, 0) - t_u_grad(0, 2);
1059 // vorticity(2) = t_u_grad(0, 1) - t_u_grad(1, 0);
1060
1061 CHKERR postProcMesh.tag_set_data(th_vorticity, &mapGaussPts[gg], 1,
1062 &vorticity(0));
1063 CHKERR postProcMesh.tag_set_data(th_q, &mapGaussPts[gg], 1, &q);
1064 CHKERR postProcMesh.tag_set_data(th_l2, &mapGaussPts[gg], 1, &l2);
1065 ++t_u_grad;
1066 }
1067
1069}
1070
1071VectorDouble3 stokes_flow_velocity(double x, double y, double z) {
1072 double r = std::sqrt(x * x + y * y + z * z);
1073 double theta = acos(x / r);
1074 double phi = atan2(y, z);
1075 double ur = cos(theta) * (1.0 + 0.5 / (r * r * r) - 1.5 / r);
1076 double ut = -sin(theta) * (1.0 - 0.25 / (r * r * r) - 0.75 / r);
1077 VectorDouble3 res(3);
1078 res[0] = ur * cos(theta) - ut * sin(theta);
1079 res[1] = ur * sin(theta) * sin(phi) + ut * cos(theta) * sin(phi);
1080 res[2] = ur * sin(theta) * cos(phi) + ut * cos(theta) * cos(phi);
1081 return res;
1082}
static Index< 'M', 3 > M
VectorDouble3 stokes_flow_velocity(double x, double y, double z)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
static double phi
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
static __CLPK_integer lapack_dsyev(char jobz, char uplo, __CLPK_integer n, __CLPK_doublereal *a, __CLPK_integer lda, __CLPK_doublereal *w, __CLPK_doublereal *work, __CLPK_integer lwork)
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
Definition Types.hpp:92
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
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.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
boost::shared_ptr< const NumeredEntFiniteElement > getNumeredEntFiniteElementPtr() const
Return raw pointer to NumeredEntFiniteElement.
Calculate inverse of jacobian for face element.
Get value at integration points for scalar field.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Operator for fat prism element updating integration weights in the volume.
Set inverse jacobian to base functions.
Transform local reference derivatives of shape functions to global derivatives.
Assemble linear (symmetric) part of the diagonal block of the LHS Operator for assembling linear (sym...
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Assemble non-linear (non-symmetric) part of the diagonal block of the LHS Operator for assembling non...
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode aSsemble(EntData &row_data, EntData &col_data)
virtual MoFEMErrorCode iNtegrate(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)
Operator for bi-linear form, usually to calculate values on left hand side.
Assemble off-diagonal block of the LHS Operator for assembling off-diagonal block of the LHS.
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
MoFEMErrorCode doWork(int row_side, EntityType row_type, EntData &row_data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode aSsemble(EntData &data)
Assemble the pressure component of the RHS vector.
MoFEMErrorCode iNtegrate(EntData &data)
Integrate local constrains vector.
Assemble linear part of the velocity component of the RHS vector.
MoFEMErrorCode iNtegrate(EntData &data)
Integrate local entity vector.
Assemble non-linear part of the velocity component of the RHS vector.
Calculate drag force on the fluid-solid interface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Calculate drag traction on the fluid-solid interface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Post processing output of drag traction on the fluid-solid interface.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
static MoFEMErrorCode setCalcDragOperators(boost::shared_ptr< FaceElementForcesAndSourcesCore > dragFe, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for calculating drag force on the solid surface.
static MoFEMErrorCode setPostProcDragOperators(boost::shared_ptr< PostProcFace > postProcDragPtr, boost::shared_ptr< VolumeElementForcesAndSourcesCoreOnSide > sideDragFe, std::string side_fe_name, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data)
Setting up operators for post processing output of drag traction.
static MoFEMErrorCode setNavierStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Navier-Stokes equations.
static MoFEMErrorCode setStokesOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > feRhs, boost::shared_ptr< VolumeElementForcesAndSourcesCore > feLhs, const std::string velocity_field, const std::string pressure_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for solving Stokes equations.
static MoFEMErrorCode setCalcVolumeFluxOperators(boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe_flux_ptr, const std::string velocity_field, boost::shared_ptr< CommonData > common_data, const EntityType type=MBTET)
Setting up operators for calculation of volume flux.