v0.13.1
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/* This file is part of MoFEM.
15 * MoFEM is free software: you can redistribute it and/or modify it under
16 * the terms of the GNU Lesser General Public License as published by the
17 * Free Software Foundation, either version 3 of the License, or (at your
18 * option) any later version.
19 *
20 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
21 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
22 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
23 * License for more details.
24 *
25 * You should have received a copy of the GNU Lesser General Public
26 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
27
28#include <MoFEM.hpp>
29using namespace MoFEM;
30using namespace boost::numeric;
33
35 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_rhs_ptr,
36 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs_ptr,
37 const std::string velocity_field, const std::string pressure_field,
38 boost::shared_ptr<CommonData> common_data, const EntityType type) {
40
41 for (auto &sit : common_data->setOfBlocksData) {
42
43 if (type == MBPRISM) {
44 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
45 fe_lhs_ptr->getOpPtrVector().push_back(
47 fe_lhs_ptr->getOpPtrVector().push_back(
48 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
49 fe_lhs_ptr->getOpPtrVector().push_back(
50 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
51 fe_rhs_ptr->getOpPtrVector().push_back(
53 fe_rhs_ptr->getOpPtrVector().push_back(
54 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
55 fe_rhs_ptr->getOpPtrVector().push_back(
56 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
57 }
58
59 fe_lhs_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
60 velocity_field, common_data->velPtr));
61 fe_lhs_ptr->getOpPtrVector().push_back(
62 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
63 common_data->gradVelPtr));
64
65 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsDiagLin(
66 velocity_field, velocity_field, common_data, sit.second));
67 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsDiagNonLin(
68 velocity_field, velocity_field, common_data, sit.second));
69 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsOffDiag(
70 velocity_field, pressure_field, common_data, sit.second));
71
72 fe_rhs_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
73 velocity_field, common_data->velPtr));
74 fe_rhs_ptr->getOpPtrVector().push_back(
75 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
76 common_data->gradVelPtr));
77 fe_rhs_ptr->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
78 pressure_field, common_data->pressPtr));
79
80 fe_rhs_ptr->getOpPtrVector().push_back(
81 new OpAssembleRhsVelocityLin(velocity_field, common_data, sit.second));
82 fe_rhs_ptr->getOpPtrVector().push_back(new OpAssembleRhsVelocityNonLin(
83 velocity_field, common_data, sit.second));
84 fe_rhs_ptr->getOpPtrVector().push_back(
85 new OpAssembleRhsPressure(pressure_field, common_data, sit.second));
86 }
87
89};
90
92 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_rhs_ptr,
93 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_lhs_ptr,
94 const std::string velocity_field, const std::string pressure_field,
95 boost::shared_ptr<CommonData> common_data, const EntityType type) {
97
98 for (auto &sit : common_data->setOfBlocksData) {
99
100 if (type == MBPRISM) {
101 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
102 fe_lhs_ptr->getOpPtrVector().push_back(
104 fe_lhs_ptr->getOpPtrVector().push_back(
105 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
106 fe_lhs_ptr->getOpPtrVector().push_back(
107 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
108 fe_rhs_ptr->getOpPtrVector().push_back(
110 fe_rhs_ptr->getOpPtrVector().push_back(
111 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
112 fe_rhs_ptr->getOpPtrVector().push_back(
113 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
114 }
115
116 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsDiagLin(
117 velocity_field, velocity_field, common_data, sit.second));
118 fe_lhs_ptr->getOpPtrVector().push_back(new OpAssembleLhsOffDiag(
119 velocity_field, pressure_field, common_data, sit.second));
120
121 fe_rhs_ptr->getOpPtrVector().push_back(
122 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
123 common_data->gradVelPtr));
124
125 fe_rhs_ptr->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
126 pressure_field, common_data->pressPtr));
127 fe_rhs_ptr->getOpPtrVector().push_back(
128 new OpAssembleRhsVelocityLin(velocity_field, common_data, sit.second));
129 fe_rhs_ptr->getOpPtrVector().push_back(
130 new OpAssembleRhsPressure(pressure_field, common_data, sit.second));
131 }
132
134};
135
137 boost::shared_ptr<VolumeElementForcesAndSourcesCore> fe_flux_ptr,
138 const std::string velocity_field, boost::shared_ptr<CommonData> common_data,
139 const EntityType type) {
141
142 for (auto &sit : common_data->setOfBlocksData) {
143
144 if (type == MBPRISM) {
145 boost::shared_ptr<MatrixDouble> inv_jac_ptr(new MatrixDouble);
146 fe_flux_ptr->getOpPtrVector().push_back(
148 fe_flux_ptr->getOpPtrVector().push_back(
149 new OpCalculateInvJacForFatPrism(inv_jac_ptr));
150 fe_flux_ptr->getOpPtrVector().push_back(
151 new OpSetInvJacH1ForFatPrism(inv_jac_ptr));
152 }
153 fe_flux_ptr->getOpPtrVector().push_back(new OpCalculateVectorFieldValues<3>(
154 velocity_field, common_data->velPtr));
155 fe_flux_ptr->getOpPtrVector().push_back(
156 new OpCalcVolumeFlux(velocity_field, common_data, sit.second));
157 }
158
160};
161
163 boost::shared_ptr<FaceElementForcesAndSourcesCore> dragFe,
164 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideDragFe,
165 std::string side_fe_name, const std::string velocity_field,
166 const std::string pressure_field,
167 boost::shared_ptr<CommonData> common_data) {
169
170 auto det_ptr = boost::make_shared<VectorDouble>();
171 auto jac_ptr = boost::make_shared<MatrixDouble>();
172 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
173
174 for (auto &sit : common_data->setOfFacesData) {
175 sideDragFe->getOpPtrVector().push_back(
176 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
177 common_data->gradVelPtr));
178 dragFe->getOpPtrVector().push_back(new OpCalculateHOJac<2>(jac_ptr));
179 dragFe->getOpPtrVector().push_back(
180 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
181 dragFe->getOpPtrVector().push_back(
182 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
183
184 dragFe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
185 pressure_field, common_data->pressPtr));
186
187 dragFe->getOpPtrVector().push_back(
189 velocity_field, sideDragFe, side_fe_name, common_data, sit.second));
190 dragFe->getOpPtrVector().push_back(new NavierStokesElement::OpCalcDragForce(
191 velocity_field, common_data, sit.second));
192 }
194};
195
197 boost::shared_ptr<PostProcFaceOnRefinedMesh> postProcDragPtr,
198 boost::shared_ptr<VolumeElementForcesAndSourcesCoreOnSide> sideDragFe,
199 std::string side_fe_name, const std::string velocity_field,
200 const std::string pressure_field,
201 boost::shared_ptr<CommonData> common_data) {
203
204 auto det_ptr = boost::make_shared<VectorDouble>();
205 auto jac_ptr = boost::make_shared<MatrixDouble>();
206 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
207
208 for (auto &sit : common_data->setOfFacesData) {
209 sideDragFe->getOpPtrVector().push_back(
210 new OpCalculateVectorFieldGradient<3, 3>(velocity_field,
211 common_data->gradVelPtr));
212 postProcDragPtr->getOpPtrVector().push_back(
213 new OpCalculateHOJac<2>(jac_ptr));
214 postProcDragPtr->getOpPtrVector().push_back(
215 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
216 postProcDragPtr->getOpPtrVector().push_back(
217 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
218
219 postProcDragPtr->getOpPtrVector().push_back(
220 new OpCalculateVectorFieldValues<3>(velocity_field,
221 common_data->velPtr));
222 postProcDragPtr->getOpPtrVector().push_back(
223 new OpCalculateScalarFieldValues(pressure_field,
224 common_data->pressPtr));
225
226 postProcDragPtr->getOpPtrVector().push_back(
228 velocity_field, sideDragFe, side_fe_name, common_data, sit.second));
229 postProcDragPtr->getOpPtrVector().push_back(
231 velocity_field, postProcDragPtr->postProcMesh,
232 postProcDragPtr->mapGaussPts, common_data, sit.second));
233 }
235};
236
238 int row_side, int col_side, EntityType row_type, EntityType col_type,
239 EntData &row_data, EntData &col_data) {
240
242
243 row_nb_dofs = row_data.getIndices().size();
244 if (!row_nb_dofs)
246 col_nb_dofs = col_data.getIndices().size();
247 if (!col_nb_dofs)
249
250 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
251 blockData.eNts.end()) {
253 }
254
255 row_nb_gauss_pts = row_data.getN().size1();
256 if (!row_nb_gauss_pts) {
258 }
259
260 isOnDiagonal = (row_type == col_type) && (row_side == col_side);
261
262 // Set size can clear local tangent matrix
263 locMat.resize(row_nb_dofs, col_nb_dofs, false);
264 locMat.clear();
265
266 CHKERR iNtegrate(row_data, col_data);
267
268 CHKERR aSsemble(row_data, col_data);
269
271}
272
274 EntData &col_data) {
276
277 const int *row_ind = &*row_data.getIndices().begin();
278 const int *col_ind = &*col_data.getIndices().begin();
279
280 Mat B = getFEMethod()->ksp_B != PETSC_NULL ? getFEMethod()->ksp_B
281 : getFEMethod()->snes_B;
282
283 CHKERR MatSetValues(B, row_nb_dofs, row_ind, col_nb_dofs, col_ind,
284 &*locMat.data().begin(), ADD_VALUES);
285
286 if (!diagonalBlock || (sYmm && !isOnDiagonal)) {
287 locMat = trans(locMat);
288 CHKERR MatSetValues(B, col_nb_dofs, col_ind, row_nb_dofs, row_ind,
289 &*locMat.data().begin(), ADD_VALUES);
290 }
291
293}
294
297 EntData &col_data) {
299
300 const int row_nb_base_functions = row_data.getN().size2();
301 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
302
306
307 // INTEGRATION
308 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
309
310 // Get volume and integration weight
311 double w = getVolume() * getGaussPts()(3, gg);
312
313 int row_bb = 0;
314 for (; row_bb != row_nb_dofs / 3; row_bb++) {
315
316 t1(i) = w * row_diff_base_functions(i);
317
318 auto base_functions = col_data.getFTensor0N(gg, 0);
319 for (int col_bb = 0; col_bb != col_nb_dofs; col_bb++) {
320
321 FTensor::Tensor1<double *, 3> k(&locMat(3 * row_bb + 0, col_bb),
322 &locMat(3 * row_bb + 1, col_bb),
323 &locMat(3 * row_bb + 2, col_bb));
324
325 k(i) -= t1(i) * base_functions;
326 ++base_functions;
327 }
328 ++row_diff_base_functions;
329 }
330 for (; row_bb != row_nb_base_functions; row_bb++) {
331 ++row_diff_base_functions;
332 }
333 }
334
336}
337
340 EntData &col_data) {
342
343 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
345 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
346 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
347 &m(r + 2, c + 2));
348 };
349
350 const int row_nb_base_functions = row_data.getN().size2();
351
352 auto row_diff_base_functions = row_data.getFTensor1DiffN<3>();
353
356
357 // integrate local matrix for entity block
358 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
359
360 // Get volume and integration weight
361 double w = getVolume() * getGaussPts()(3, gg);
362 double const alpha = w * blockData.viscousCoef;
363
364 int row_bb = 0;
365 for (; row_bb != row_nb_dofs / 3; row_bb++) {
366
367 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
368
369 const int final_bb = isOnDiagonal ? row_bb + 1 : col_nb_dofs / 3;
370
371 int col_bb = 0;
372
373 for (; col_bb != final_bb; col_bb++) {
374
375 auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
376
377 for (int i : {0, 1, 2}) {
378 for (int j : {0, 1, 2}) {
379 t_assemble(i, i) +=
380 alpha * row_diff_base_functions(j) * col_diff_base_functions(j);
381 t_assemble(i, j) +=
382 alpha * row_diff_base_functions(j) * col_diff_base_functions(i);
383 }
384 }
385
386 // Next base function for column
387 ++col_diff_base_functions;
388 }
389
390 // Next base function for row
391 ++row_diff_base_functions;
392 }
393
394 for (; row_bb != row_nb_base_functions; row_bb++) {
395 ++row_diff_base_functions;
396 }
397 }
398
399 if (isOnDiagonal) {
400 for (int row_bb = 0; row_bb != row_nb_dofs / 3; row_bb++) {
401 int col_bb = 0;
402 for (; col_bb != row_bb + 1; col_bb++) {
403 auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
404 auto t_off_side = get_tensor2(locMat, 3 * col_bb, 3 * row_bb);
405 t_off_side(i, j) = t_assemble(j, i);
406 }
407 }
408 }
409
411}
412
415 EntData &col_data) {
417
418 auto get_tensor2 = [](MatrixDouble &m, const int r, const int c) {
420 &m(r + 0, c + 0), &m(r + 0, c + 1), &m(r + 0, c + 2), &m(r + 1, c + 0),
421 &m(r + 1, c + 1), &m(r + 1, c + 2), &m(r + 2, c + 0), &m(r + 2, c + 1),
422 &m(r + 2, c + 2));
423 };
424
425 const int row_nb_base_functions = row_data.getN().size2();
426
427 auto row_base_functions = row_data.getFTensor0N();
428
429 auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
430 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
431
434
435 // integrate local matrix for entity block
436 for (int gg = 0; gg != row_nb_gauss_pts; gg++) {
437
438 // Get volume and integration weight
439 double w = getVolume() * getGaussPts()(3, gg);
440 double const beta = w * blockData.inertiaCoef;
441
442 int row_bb = 0;
443 for (; row_bb != row_nb_dofs / 3; row_bb++) {
444
445 auto col_diff_base_functions = col_data.getFTensor1DiffN<3>(gg, 0);
446 auto col_base_functions = col_data.getFTensor0N(gg, 0);
447
448 const int final_bb = col_nb_dofs / 3;
449 int col_bb = 0;
450
451 for (; col_bb != final_bb; col_bb++) {
452
453 auto t_assemble = get_tensor2(locMat, 3 * row_bb, 3 * col_bb);
454
455 for (int i : {0, 1, 2}) {
456 for (int j : {0, 1, 2}) {
457 t_assemble(i, j) +=
458 beta * col_base_functions * t_u_grad(i, j) * row_base_functions;
459 t_assemble(i, i) +=
460 beta * t_u(j) * col_diff_base_functions(j) * row_base_functions;
461 }
462 }
463
464 // Next base function for column
465 ++col_diff_base_functions;
466 ++col_base_functions;
467 }
468
469 // Next base function for row
470 ++row_base_functions;
471 }
472
473 for (; row_bb != row_nb_base_functions; row_bb++) {
474 ++row_base_functions;
475 }
476
477 ++t_u;
478 ++t_u_grad;
479 }
480
482}
483
485 EntityType row_type,
486 EntData &row_data) {
488 // get number of dofs on row
489 nbRows = row_data.getIndices().size();
490 if (!nbRows)
492 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
493 blockData.eNts.end()) {
495 }
496 // get number of integration points
497 nbIntegrationPts = getGaussPts().size2();
498
499 // integrate local vector
500 CHKERR iNtegrate(row_data);
501 // assemble local vector
502 CHKERR aSsemble(row_data);
504}
505
508 // get global indices of local vector
509 const int *indices = &*data.getIndices().data().begin();
510 // get values from local vector
511 const double *vals = &*locVec.data().begin();
512 Vec f = getFEMethod()->ksp_f != PETSC_NULL ? getFEMethod()->ksp_f
513 : getFEMethod()->snes_f;
514 // assemble vector
515 CHKERR VecSetValues(f, nbRows, indices, vals, ADD_VALUES);
516
518}
519
523
524 auto get_tensor1 = [](VectorDouble &v, const int r) {
526 &v(r + 0), &v(r + 1), &v(r + 2));
527 };
528
529 // set size of local vector
530 locVec.resize(nbRows, false);
531 // clear local entity vector
532 locVec.clear();
533
534 int nb_base_functions = data.getN().size2();
535
536 // get base function gradient on rows
537 auto t_v_grad = data.getFTensor1DiffN<3>();
538
539 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
540 auto t_p = getFTensor0FromVec(*commonData->pressPtr);
541
544
545 // loop over all integration points
546 for (int gg = 0; gg != nbIntegrationPts; gg++) {
547
548 double w = getVolume() * getGaussPts()(3, gg);
549
550 // evaluate constant term
551 const double alpha = w * blockData.viscousCoef;
552
553 auto t_a = get_tensor1(locVec, 0);
554 int rr = 0;
555
556 // loop over base functions
557 for (; rr != nbRows / 3; rr++) {
558 // add to local vector source term
559
560 t_a(i) += alpha * t_u_grad(i, j) * t_v_grad(j);
561 t_a(j) += alpha * t_u_grad(i, j) * t_v_grad(i);
562
563 t_a(i) -= w * t_p * t_v_grad(i);
564
565 ++t_a; // move to next element of local vector
566 ++t_v_grad; // move to next gradient of base function
567 }
568
569 for (; rr != nb_base_functions; rr++) {
570 ++t_v_grad;
571 }
572
573 ++t_u_grad;
574 ++t_p;
575 }
576
578}
579
583
584 auto get_tensor1 = [](VectorDouble &v, const int r) {
586 &v(r + 0), &v(r + 1), &v(r + 2));
587 };
588
589 // set size of local vector
590 locVec.resize(nbRows, false);
591 // clear local entity vector
592 locVec.clear();
593
594 // get base functions on entity
595 auto t_v = data.getFTensor0N();
596
597 int nb_base_functions = data.getN().size2();
598
599 auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
600 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
601
604
605 // loop over all integration points
606 for (int gg = 0; gg != nbIntegrationPts; gg++) {
607
608 double w = getVolume() * getGaussPts()(3, gg);
609 // evaluate constant term
610 const double beta = w * blockData.inertiaCoef;
611
612 auto t_a = get_tensor1(locVec, 0);
613 int rr = 0;
614
615 // loop over base functions
616 for (; rr != nbRows / 3; rr++) {
617 // add to local vector source term
618
619 t_a(i) += beta * t_v * t_u_grad(i, j) * t_u(j);
620
621 ++t_a; // move to next element of local vector
622 ++t_v; // move to next base function
623 }
624
625 for (; rr != nb_base_functions; rr++) {
626 ++t_v;
627 }
628
629 ++t_u;
630 ++t_u_grad;
631 }
633}
634
638
639 // commonData->getBlockData(blockData);
640
641 // set size to local vector
642 locVec.resize(nbRows, false);
643 // clear local vector
644 locVec.clear();
645 // get base function
646 auto t_q = data.getFTensor0N();
647
648 int nb_base_functions = data.getN().size2();
649 // get solution at integration point
650
651 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
652
653 // FTensor::Index<'i', 3> i;
654
655 // make loop over integration points
656 for (int gg = 0; gg != nbIntegrationPts; gg++) {
657 // evaluate function on boundary and scale it by area and integration
658 // weight
659
660 double w = getVolume() * getGaussPts()(3, gg);
661
662 // get element of vector
664 int rr = 0;
665 for (; rr != nbRows; rr++) {
666
667 for (int ii : {0, 1, 2}) {
668 t_a -= w * t_q * t_u_grad(ii, ii);
669 }
670
671 ++t_a;
672 ++t_q;
673 }
674
675 for (; rr != nb_base_functions; rr++) {
676 ++t_q;
677 }
678
679 ++t_u_grad;
680 }
682}
683
686 EntData &data) {
688
689 if (type != MBVERTEX)
690 PetscFunctionReturn(0);
691
692 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
693 blockData.eNts.end()) {
695 }
696
697 const int nb_gauss_pts = getGaussPts().size2();
698
699 auto t_u = getFTensor1FromMat<3>(*commonData->velPtr);
700 auto t_w = getFTensor0IntegrationWeight(); ///< Integration weight
701
703
705 t_flux(i) = 0.0; // Zero entries
706
707 // loop over all integration points
708 for (int gg = 0; gg != nb_gauss_pts; ++gg) {
709
710 double vol = getVolume();
711 t_flux(i) += t_w * vol * t_u(i);
712
713 ++t_w;
714 ++t_u;
715 }
716
717 // Set array of indices
718 constexpr std::array<int, 3> indices = {0, 1, 2};
719
720 // Assemble volumetric flux
721 CHKERR VecSetValues(commonData->volumeFluxVec, 3, indices.data(), &t_flux(0),
722 ADD_VALUES);
723
725}
726
729 EntData &data) {
731 if (type != MBVERTEX)
732 PetscFunctionReturn(0);
733
734 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
735 blockData.eNts.end()) {
737 }
738
739 const int nb_gauss_pts = getGaussPts().size2();
740
741 auto pressure_drag_at_gauss_pts =
742 getFTensor1FromMat<3>(*commonData->pressureDragTract);
743 auto shear_drag_at_gauss_pts =
744 getFTensor1FromMat<3>(*commonData->shearDragTract);
745 auto total_drag_at_gauss_pts =
746 getFTensor1FromMat<3>(*commonData->totalDragTract);
747
748 FTensor::Tensor1<double, 3> t_pressure_drag_force;
749 FTensor::Tensor1<double, 3> t_shear_drag_force;
750 FTensor::Tensor1<double, 3> t_total_drag_force;
751
753
754 t_pressure_drag_force(i) = 0.0; // Zero entries
755 t_shear_drag_force(i) = 0.0; // Zero entries
756 t_total_drag_force(i) = 0.0; // Zero entries
757
758 auto t_w = getFTensor0IntegrationWeight();
759 auto t_normal = getFTensor1NormalsAtGaussPts();
760
761 for (int gg = 0; gg != nb_gauss_pts; gg++) {
762
763 double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
764 double w = t_w * nrm2 * 0.5;
765
766 t_pressure_drag_force(i) += w * pressure_drag_at_gauss_pts(i);
767 t_shear_drag_force(i) += w * shear_drag_at_gauss_pts(i);
768 t_total_drag_force(i) += w * total_drag_at_gauss_pts(i);
769
770 ++t_w;
771 ++t_normal;
772 ++pressure_drag_at_gauss_pts;
773 ++shear_drag_at_gauss_pts;
774 ++total_drag_at_gauss_pts;
775 }
776
777 // Set array of indices
778 constexpr std::array<int, 3> indices = {0, 1, 2};
779
780 CHKERR VecSetValues(commonData->pressureDragForceVec, 3, indices.data(),
781 &t_pressure_drag_force(0), ADD_VALUES);
782 CHKERR VecSetValues(commonData->shearDragForceVec, 3, indices.data(),
783 &t_shear_drag_force(0), ADD_VALUES);
784 CHKERR VecSetValues(commonData->totalDragForceVec, 3, indices.data(),
785 &t_total_drag_force(0), ADD_VALUES);
786
788}
789
792 EntData &data) {
794
795 if (type != MBVERTEX)
796 PetscFunctionReturn(0);
797
798 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
799 blockData.eNts.end()) {
801 }
802
803 CHKERR loopSideVolumes(sideFeName, *sideFe);
804
805 const int nb_gauss_pts = getGaussPts().size2();
806
807 auto t_p = getFTensor0FromVec(*commonData->pressPtr);
808 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
809
810 auto t_normal = getFTensor1NormalsAtGaussPts();
811
812 FTensor::Tensor1<double, 3> t_unit_normal;
813
816
817 commonData->pressureDragTract->resize(3, nb_gauss_pts, false);
818 commonData->pressureDragTract->clear();
819
820 commonData->shearDragTract->resize(3, nb_gauss_pts, false);
821 commonData->shearDragTract->clear();
822
823 commonData->totalDragTract->resize(3, nb_gauss_pts, false);
824 commonData->totalDragTract->clear();
825
826 auto pressure_drag_at_gauss_pts =
827 getFTensor1FromMat<3>(*commonData->pressureDragTract);
828 auto shear_drag_at_gauss_pts =
829 getFTensor1FromMat<3>(*commonData->shearDragTract);
830 auto total_drag_at_gauss_pts =
831 getFTensor1FromMat<3>(*commonData->totalDragTract);
832
833 for (int gg = 0; gg != nb_gauss_pts; gg++) {
834
835 double nrm2 = std::sqrt(t_normal(i) * t_normal(i));
836 t_unit_normal(i) = t_normal(i) / nrm2;
837
838 double mu = blockData.fluidViscosity;
839
840 pressure_drag_at_gauss_pts(i) = t_p * t_unit_normal(i);
841 shear_drag_at_gauss_pts(i) =
842 -mu * (t_u_grad(i, j) + t_u_grad(j, i)) * t_unit_normal(j);
843 total_drag_at_gauss_pts(i) =
844 pressure_drag_at_gauss_pts(i) + shear_drag_at_gauss_pts(i);
845
846 ++pressure_drag_at_gauss_pts;
847 ++shear_drag_at_gauss_pts;
848 ++total_drag_at_gauss_pts;
849 ++t_p;
850 ++t_u_grad;
851 ++t_normal;
852 }
854}
855
858 EntData &data) {
860 if (type != MBVERTEX)
861 PetscFunctionReturn(0);
862
863 double def_VAL[9];
864 bzero(def_VAL, 9 * sizeof(double));
865
866 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
867 blockData.eNts.end()) {
869 }
870
871 Tag th_velocity;
872 Tag th_pressure;
873 Tag th_velocity_grad;
874 Tag th_shear_drag;
875 Tag th_pressure_drag;
876 Tag th_total_drag;
877
878 CHKERR postProcMesh.tag_get_handle("VELOCITY", 3, MB_TYPE_DOUBLE, th_velocity,
879 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
880 CHKERR postProcMesh.tag_get_handle("PRESSURE", 1, MB_TYPE_DOUBLE, th_pressure,
881 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
882 CHKERR postProcMesh.tag_get_handle("VELOCITY_GRAD", 9, MB_TYPE_DOUBLE,
883 th_velocity_grad,
884 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
885
886 CHKERR postProcMesh.tag_get_handle("PRESSURE_DRAG", 3, MB_TYPE_DOUBLE,
887 th_pressure_drag,
888 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
889 CHKERR postProcMesh.tag_get_handle("SHEAR_DRAG", 3, MB_TYPE_DOUBLE,
890 th_shear_drag,
891 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
892 CHKERR postProcMesh.tag_get_handle("TOTAL_DRAG", 3, MB_TYPE_DOUBLE,
893 th_total_drag,
894 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
895
896 auto t_p = getFTensor0FromVec(*commonData->pressPtr);
897
898 const int nb_gauss_pts = commonData->pressureDragTract->size2();
899
900 MatrixDouble velGradMat;
901 velGradMat.resize(3, 3);
902 VectorDouble velVec;
903 velVec.resize(3);
904 VectorDouble pressDragVec;
905 pressDragVec.resize(3);
906 VectorDouble viscDragVec;
907 viscDragVec.resize(3);
908 VectorDouble totDragVec;
909 totDragVec.resize(3);
910
911 for (int gg = 0; gg != nb_gauss_pts; gg++) {
912
913 for (int i : {0, 1, 2}) {
914 for (int j : {0, 1, 2}) {
915 velGradMat(i, j) = (*commonData->gradVelPtr)(i * 3 + j, gg);
916 }
917 velVec(i) = (*commonData->velPtr)(i, gg);
918 pressDragVec(i) = (*commonData->pressureDragTract)(i, gg);
919 viscDragVec(i) = (*commonData->shearDragTract)(i, gg);
920 totDragVec(i) = (*commonData->totalDragTract)(i, gg);
921 }
922 CHKERR postProcMesh.tag_set_data(th_velocity, &mapGaussPts[gg], 1,
923 &velVec(0));
924 CHKERR postProcMesh.tag_set_data(th_pressure, &mapGaussPts[gg], 1, &t_p);
925 CHKERR postProcMesh.tag_set_data(th_velocity_grad, &mapGaussPts[gg], 1,
926 &velGradMat(0, 0));
927
928 CHKERR postProcMesh.tag_set_data(th_pressure_drag, &mapGaussPts[gg], 1,
929 &pressDragVec(0));
930
931 CHKERR postProcMesh.tag_set_data(th_shear_drag, &mapGaussPts[gg], 1,
932 &viscDragVec(0));
933
934 CHKERR postProcMesh.tag_set_data(th_total_drag, &mapGaussPts[gg], 1,
935 &totDragVec(0));
936
937 ++t_p;
938 }
939
941}
942
945 EntData &data) {
947 if (type != MBVERTEX)
948 PetscFunctionReturn(0);
949 double def_VAL[9];
950 bzero(def_VAL, 9 * sizeof(double));
951
952 if (blockData.eNts.find(getNumeredEntFiniteElementPtr()->getEnt()) ==
953 blockData.eNts.end()) {
955 }
956 // commonData->getBlockData(blockData);
957
958 Tag th_vorticity;
959 Tag th_q;
960 Tag th_l2;
961 CHKERR postProcMesh.tag_get_handle("VORTICITY", 3, MB_TYPE_DOUBLE,
962 th_vorticity, MB_TAG_CREAT | MB_TAG_SPARSE,
963 def_VAL);
964 CHKERR postProcMesh.tag_get_handle("Q", 1, MB_TYPE_DOUBLE, th_q,
965 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
966 CHKERR postProcMesh.tag_get_handle("L2", 1, MB_TYPE_DOUBLE, th_l2,
967 MB_TAG_CREAT | MB_TAG_SPARSE, def_VAL);
968
969 auto t_u_grad = getFTensor2FromMat<3, 3>(*commonData->gradVelPtr);
970 // auto p = getFTensor0FromVec(*commonData->pressPtr);
971
972 const int nb_gauss_pts = commonData->gradVelPtr->size2();
973 // const int nb_gauss_pts2 = commonData->pressPtr->size();
974
975 // const double lambda = commonData->lAmbda;
976 // const double mu = commonData->mU;
977
978 // FTensor::Index<'i', 3> i;
979 // FTensor::Index<'j', 3> j;
980 // FTensor::Index<'j', 3> k;
982 // FTensor::Tensor2<double,3,3> t_s;
983 double q;
984 double l2;
985 // FTensor::Tensor2<double, 3, 3> stress;
986 MatrixDouble S;
987 MatrixDouble Omega;
989
990 S.resize(3, 3);
991 Omega.resize(3, 3);
992 M.resize(3, 3);
993
994 for (int gg = 0; gg != nb_gauss_pts; gg++) {
995
996 vorticity(0) = t_u_grad(2, 1) - t_u_grad(1, 2);
997 vorticity(1) = t_u_grad(0, 2) - t_u_grad(2, 0);
998 vorticity(2) = t_u_grad(1, 0) - t_u_grad(0, 1);
999
1000 q = 0;
1001 for (int i = 0; i != 3; i++) {
1002 for (int j = 0; j != 3; j++) {
1003 q += -0.5 * t_u_grad(i, j) * t_u_grad(j, i);
1004 }
1005 }
1006 for (int i = 0; i != 3; i++) {
1007 for (int j = 0; j != 3; j++) {
1008 S(i, j) = 0.5 * (t_u_grad(i, j) + t_u_grad(j, i));
1009 Omega(i, j) = 0.5 * (t_u_grad(i, j) - t_u_grad(j, i));
1010 M(i, j) = 0.0;
1011 }
1012 }
1013
1014 for (int i = 0; i != 3; i++) {
1015 for (int j = 0; j != 3; j++) {
1016 for (int k = 0; k != 3; k++) {
1017 M(i, j) += S(i, k) * S(k, j) + Omega(i, k) * Omega(k, j);
1018 }
1019 }
1020 }
1021
1022 MatrixDouble eigen_vectors = M;
1023 VectorDouble eigen_values(3);
1024
1025 // LAPACK - eigenvalues and vectors. Applied twice for initial creates
1026 // memory space
1027 int n = 3, lda = 3, info, lwork = -1;
1028 double wkopt;
1029 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
1030 &(eigen_values.data()[0]), &wkopt, lwork);
1031 if (info != 0)
1032 SETERRQ1(PETSC_COMM_SELF, 1,
1033 "is something wrong with lapack_dsyev info = %d", info);
1034 lwork = (int)wkopt;
1035 double work[lwork];
1036 info = lapack_dsyev('V', 'U', n, &(eigen_vectors.data()[0]), lda,
1037 &(eigen_values.data()[0]), work, lwork);
1038 if (info != 0)
1039 SETERRQ1(PETSC_COMM_SELF, 1,
1040 "is something wrong with lapack_dsyev info = %d", info);
1041
1042 map<double, int> eigen_sort;
1043 eigen_sort[eigen_values[0]] = 0;
1044 eigen_sort[eigen_values[1]] = 1;
1045 eigen_sort[eigen_values[2]] = 2;
1046
1047 // prin_stress_vect.clear();
1048 VectorDouble prin_vals_vect(3);
1049 prin_vals_vect.clear();
1050
1051 int ii = 0;
1052 for (map<double, int>::reverse_iterator mit = eigen_sort.rbegin();
1053 mit != eigen_sort.rend(); mit++) {
1054 prin_vals_vect[ii] = eigen_values[mit->second];
1055 // for (int dd = 0; dd != 3; dd++) {
1056 // prin_stress_vect(ii, dd) = eigen_vectors.data()[3 * mit->second +
1057 // dd];
1058 // }
1059 ii++;
1060 }
1061
1062 l2 = prin_vals_vect[1];
1063 // cout << prin_vals_vect << endl;
1064 // cout << "-0.5 sum: " << -0.5 * (prin_vals_vect[0] + prin_vals_vect[1]
1065 // + prin_vals_vect[2]) << endl; cout << "q: " << q << endl;
1066
1067 // t_s(i,j) = 0.5*(t)
1068
1069 // vorticity(0) = t_u_grad(1, 2) - t_u_grad(2, 1);
1070 // vorticity(1) = t_u_grad(2, 0) - t_u_grad(0, 2);
1071 // vorticity(2) = t_u_grad(0, 1) - t_u_grad(1, 0);
1072
1073 CHKERR postProcMesh.tag_set_data(th_vorticity, &mapGaussPts[gg], 1,
1074 &vorticity(0));
1075 CHKERR postProcMesh.tag_set_data(th_q, &mapGaussPts[gg], 1, &q);
1076 CHKERR postProcMesh.tag_set_data(th_l2, &mapGaussPts[gg], 1, &l2);
1077 ++t_u_grad;
1078 }
1079
1081}
1082
1083VectorDouble3 stokes_flow_velocity(double x, double y, double z) {
1084 double r = std::sqrt(x * x + y * y + z * z);
1085 double theta = acos(x / r);
1086 double phi = atan2(y, z);
1087 double ur = cos(theta) * (1.0 + 0.5 / (r * r * r) - 1.5 / r);
1088 double ut = -sin(theta) * (1.0 - 0.25 / (r * r * r) - 0.75 / r);
1089 VectorDouble3 res(3);
1090 res[0] = ur * cos(theta) - ut * sin(theta);
1091 res[1] = ur * sin(theta) * sin(phi) + ut * cos(theta) * sin(phi);
1092 res[2] = ur * sin(theta) * cos(phi) + ut * cos(theta) * cos(phi);
1093 return res;
1094}
static Index< 'M', 3 > M
static Index< 'n', 3 > n
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()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
constexpr double beta
constexpr double alpha
FTensor::Index< 'i', SPACE_DIM > i
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)
Definition: lapack_wrap.h:273
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:103
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:149
double w(const double g, const double t)
const double r
rate factor
static double phi
constexpr double mu
FTensor::Index< 'm', 3 > m
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.
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< PostProcFaceOnRefinedMesh > 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.