v0.14.0
Loading...
Searching...
No Matches
fluid_structure_eigenproblem.cpp
Go to the documentation of this file.
1/**
2 * \file eigen_elastic.cpp
3 * \FSI eigen_elastic.cpp
4 *
5 * Calculate natural frequencies in 2d and 3d problems.
6 *
7 */
8
9#include <MoFEM.hpp>
10#undef EPS
11#include <slepceps.h>
12#include <slepcpep.h>
13
14using namespace MoFEM;
15
16template <int DIM> struct ElementsAndOps {};
17
18template <> struct ElementsAndOps<2> {
22
25
28};
29
30template <> struct ElementsAndOps<3> {
34
37
40};
41
42constexpr int SPACE_DIM =
43 EXECUTABLE_DIMENSION; //< Space dimension of problem, mesh
44
49
52
55
57 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
60
67
68static double penalty = 1e2;
69static double phi = -1;
70
71// Length unit are milometer
72// Time unit miliseconds
73// Force unit Newton
74
75double rho_solid = 2.28;
76double K_solid = 95 * 1e3;
77double G_solid = 63 * 1e3;
78
79double rho_fluid = 0.998;
80double K_fluid = 2.2 * 1e3;
81double G_fluid = 1e-6;
82double MU_fluid = 1e-3;
83
84int order = 1;
85
86double eig_threshold = 1e-1;
87
88// FIXME
89auto poisson_ratio = 0.265;
90#include <OpPostProcElastic.hpp>
91
92using namespace Tutorial;
93
94/**
95 * @brief Operator tp collect data from elements on the side of Edge/Face
96 *
97 */
99
100 OpCalculateSideData(std::string field_name, std::string col_field_name,
101 UserDataOperator::OpType type);
102
103 MoFEMErrorCode doWork(int side, EntityType type, EntData &data);
104};
105
106/**
107 * @brief Operator the left hand side matrix
108 *
109 */
111public:
112 /**
113 * @brief Construct a new OpH1LhsSkeleton
114 *
115 * @param side_fe_ptr pointer to FE to evaluate side elements
116 */
117 OpLhsSkeleton(boost::shared_ptr<FaceSideEle> side_fe_ptr,
118 boost::shared_ptr<MatrixDouble> d_mat_fluid_ptr,
119 boost::shared_ptr<MatrixDouble> d_mat_solid_ptr,
120 boost::shared_ptr<MatrixDouble> d_mat_visc_fluid_ptr,
121 boost::shared_ptr<MatrixDouble> d_mat_visc_solid_ptr,
122 boost::shared_ptr<Range> r_visc_fluid_ptr,
123 boost::shared_ptr<Range> r_visc_solid_ptr, SmartPetscObj<Mat> K,
125
126 MoFEMErrorCode doWork(int side, EntityType type,
128
129private:
130 boost::shared_ptr<FaceSideEle>
131 sideFEPtr; ///< pointer to element to get data on edge/face sides
135 boost::shared_ptr<MatrixDouble> dMatFluidPtr;
136 boost::shared_ptr<MatrixDouble> dMatSolidPtr;
137 boost::shared_ptr<MatrixDouble> dViscFluidPtr;
138 boost::shared_ptr<MatrixDouble> dViscSolidPtr;
139 boost::shared_ptr<Range> rFluidPtr;
140 boost::shared_ptr<Range> rSolidPtr;
144};
145
146struct FSI {
147
148 FSI(MoFEM::Interface &m_field) : mField(m_field) {}
149
151
152private:
154
162
163 boost::shared_ptr<MatrixDouble> matGradPtr;
164 boost::shared_ptr<MatrixDouble> matStrainPtr;
165 boost::shared_ptr<MatrixDouble> matStressPtr;
166 boost::shared_ptr<MatrixDouble> matDSolidPtr;
167 boost::shared_ptr<MatrixDouble> matDFluidPtr;
168 boost::shared_ptr<MatrixDouble> matDFViscFluidPtr;
169 boost::shared_ptr<MatrixDouble> matDFViscSolidPtr;
170
174
177
180
182};
183
184//! [Create common data]
187
188 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho_solid", &rho_solid,
189 PETSC_NULL);
190 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho_fluid", &rho_fluid,
191 PETSC_NULL);
192 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-bulk_modulus_solid", &K_solid,
193 PETSC_NULL);
194 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-sheer_modulus_solid", &G_solid,
195 PETSC_NULL);
196 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-bulk_modulus_fluid", &K_fluid,
197 PETSC_NULL);
198 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-viscosity_fluid", &MU_fluid,
199 PETSC_NULL);
200 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-penalty", &penalty,
201 PETSC_NULL);
202 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-phi", &phi, PETSC_NULL);
203
204 MOFEM_LOG("FSI", Sev::inform) << "-rho_solid " << rho_solid;
205 MOFEM_LOG("FSI", Sev::inform) << "-rho_fluid " << rho_fluid;
206 MOFEM_LOG("FSI", Sev::inform) << "-bulk_modulus_solid " << K_solid;
207 MOFEM_LOG("FSI", Sev::inform) << "-sheer_modulus_solid " << G_solid;
208 MOFEM_LOG("FSI", Sev::inform) << "-bulk_modulus_fluid " << K_fluid;
209 MOFEM_LOG("FSI", Sev::inform) << "-viscosity_fluid " << MU_fluid;
210 MOFEM_LOG("FSI", Sev::inform) << "-penalty " << penalty;
211 MOFEM_LOG("FSI", Sev::inform) << "-phi " << phi;
212
213 const char *list_solvers[]{"eps", "pep"};
215 CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-solver_type", list_solvers,
216 LASTSOLVER, &solverType, PETSC_NULL);
217
218 auto set_matrial_stiffens = [&](auto &mat_d, auto K, auto G) {
219 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
220 mat_d.resize(size_symm * size_symm, 1);
221
224 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(mat_d);
225
226 constexpr double A = 1; // plane strain
227
228 // double A = (SPACE_DIM == 2) ? 2 * G / (K + (4. / 3.) * G) : 1;
229
230 t_D(i, j, k, l) = 2 * G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
231 A * (K - (2. / 3.) * G) * t_kd(i, j) * t_kd(k, l);
232
234 };
235
236 auto set_viscosity = [&](auto &mat_d, auto MU) {
237 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
238 mat_d.resize(size_symm * size_symm, 1);
239
242 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(mat_d);
243
244 t_D(i, j, k, l) = MU * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
245
247 };
248
249 matGradPtr = boost::make_shared<MatrixDouble>();
250 matStrainPtr = boost::make_shared<MatrixDouble>();
251 matStressPtr = boost::make_shared<MatrixDouble>();
252 matDSolidPtr = boost::make_shared<MatrixDouble>();
253 matDFluidPtr = boost::make_shared<MatrixDouble>();
254 matDFViscFluidPtr = boost::make_shared<MatrixDouble>();
255 matDFViscSolidPtr = boost::make_shared<MatrixDouble>();
256
257 CHKERR set_matrial_stiffens(*matDSolidPtr, K_solid, G_solid);
258 CHKERR set_matrial_stiffens(*matDFluidPtr, K_fluid, G_fluid);
259 // CHKERR set_matrial_stiffens(*matDFViscFluidPtr, 0, MU_fluid);
260 CHKERR set_viscosity(*matDFViscFluidPtr, MU_fluid);
261
262 // No solid damping
263 matDFViscSolidPtr->resize(matDFViscFluidPtr->size1(),
264 matDFViscFluidPtr->size2());
265 matDFViscSolidPtr->clear();
266
268}
269//! [Create common data]
270
271//! [Run problem]
283//! [Run problem]
284
285//! [Read mesh]
289
290 MOFEM_LOG("FSI", Sev::inform)
291 << "Read mesh for problem in " << EXECUTABLE_DIMENSION;
292 CHKERR simple->getOptions();
293
294 simple->getAddSkeletonFE() = true;
295 CHKERR simple->loadFile();
296
298
299 const std::string block_names[] = {"SOLID", "FLUID", "INTERFACE"};
300
301 if (it->getName().compare(0, block_names[0].length(), block_names[0]) ==
302 0) {
303 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, rSolid,
304 true);
305 rSolid = rSolid.subset_by_dimension(SPACE_DIM);
306 } else if (it->getName().compare(0, block_names[1].length(),
307 block_names[1]) == 0) {
308 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, rFluid,
309 true);
310 rFluid = rFluid.subset_by_dimension(SPACE_DIM);
311 } else if (it->getName().compare(0, block_names[2].length(),
312 block_names[2]) == 0) {
313 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, rInterface,
314 true);
315 rInterface = rInterface.subset_by_dimension(SPACE_DIM - 1);
316 }
317 }
318
319 for (auto d = SPACE_DIM - 1; d >= 0; --d) {
320 auto add_adj = [&](auto &r, auto dim) {
322 Range adj;
323 CHKERR mField.get_moab().get_adjacencies(
324 r.subset_by_dimension(dim), d, false, adj, moab::Interface::UNION);
325 r.merge(adj);
327 };
328 CHKERR add_adj(rSolid, SPACE_DIM);
329 CHKERR add_adj(rFluid, SPACE_DIM);
330 if (d < (SPACE_DIM - 1)) {
331 CHKERR add_adj(rInterface, SPACE_DIM - 1);
332 }
333 }
334
335 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(rSolid);
336 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(rFluid);
337 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(rInterface);
338
339 MOFEM_LOG("FSI", Sev::noisy) << "Solid " << rSolid << endl;
340 MOFEM_LOG("FSI", Sev::noisy) << "Fluid " << rFluid << endl;
341 MOFEM_LOG("FSI", Sev::noisy) << "Interface " << rInterface << endl;
342
344}
345//! [Read mesh]
346
347//! [Set up problem]
349 auto *simple = mField.getInterface<Simple>();
351 // Add field
352 CHKERR simple->addDomainField("Us", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
353 CHKERR simple->addDomainField("Uf", H1, AINSWORTH_LEGENDRE_BASE, SPACE_DIM);
354 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
355 CHKERR simple->setFieldOrder("Us", order);
356 CHKERR simple->setFieldOrder("Uf", order);
357
358 if (solverType == PEP_SOLVER) {
359
360 CHKERR mField.add_finite_element("INTERFACE");
361 auto add_fe_field = [&](auto field) {
367 };
368 CHKERR add_fe_field("Uf");
369 CHKERR add_fe_field("Us");
370
372 "INTERFACE");
373 simple->getOtherFiniteElements().push_back("INTERFACE");
374 }
375
376 CHKERR simple->defineFiniteElements();
377 CHKERR simple->setSkeletonAdjacency();
378 CHKERR simple->defineProblem();
379 CHKERR simple->buildFields();
380 CHKERR simple->buildFiniteElements();
381 CHKERR simple->buildProblem();
382
383 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
384 simple->getProblemName(), "Uf", subtract(rSolid, rInterface));
385 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
386 simple->getProblemName(), "Us", subtract(rFluid, rInterface));
387
389}
390//! [Set up problem]
391
392//! [Boundary condition]
395 auto bc_mng = mField.getInterface<BcManager>();
397
398 auto bc_mng = mField.getInterface<BcManager>();
399 auto &bc_map = bc_mng->getBcMapByBlockName();
400
401 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
402 "Uf", 0, 0);
403 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
404 "Uf", 1, 1);
405 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
406 "Uf", 2, 2);
407 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
408 "Uf", 0, 3);
409
410 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_X",
411 "Us", 0, 0);
412 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Y",
413 "Us", 1, 1);
414 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_Z",
415 "Us", 2, 2);
416 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "FIX_ALL",
417 "Us", 0, 3);
418
420}
421//! [Boundary condition]
422
423//! [Push operators to pipeline]
427 auto pipeline_mng = mField.getInterface<PipelineManager>();
428 auto bc_mng = mField.getInterface<BcManager>();
429 auto &bc_map = bc_mng->getBcMapByBlockName();
430
431 auto dm = simple->getDM();
433 M = smartMatDuplicate(K, MAT_SHARE_NONZERO_PATTERN);
434 C = smartMatDuplicate(K, MAT_SHARE_NONZERO_PATTERN);
435
436 CHKERR MatZeroEntries(K);
437 CHKERR MatZeroEntries(M);
438 CHKERR MatZeroEntries(C);
439
440 auto ptr_r_solid =
441 boost::make_shared<Range>(rSolid.subset_by_dimension(SPACE_DIM));
442 auto ptr_r_fluid =
443 boost::make_shared<Range>(rFluid.subset_by_dimension(SPACE_DIM));
444
445 auto basic_ops = [&]() {
447 auto det_ptr = boost::make_shared<VectorDouble>();
448 auto jac_ptr = boost::make_shared<MatrixDouble>();
449 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
450 pipeline_mng->getOpDomainLhsPipeline().push_back(
451 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
452 pipeline_mng->getOpDomainLhsPipeline().push_back(
453 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
454 pipeline_mng->getOpDomainLhsPipeline().push_back(
455 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
456 pipeline_mng->getOpDomainLhsPipeline().push_back(
459 };
460
461 auto integration_rule = [](int, int, int approx_order) {
462 return 2 * approx_order + 1;
463 };
464
465 auto calculate_stiffness_matrix = [&]() {
467 pipeline_mng->getDomainLhsFE().reset();
468 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
469
470 CHKERR basic_ops();
471
472 pipeline_mng->getOpDomainLhsPipeline().push_back(
473 new OpK("Us", "Us", matDSolidPtr, ptr_r_solid));
474 pipeline_mng->getOpDomainLhsPipeline().push_back(
475 new OpK("Uf", "Uf", matDFluidPtr, ptr_r_fluid));
476
477 pipeline_mng->getDomainLhsFE()->B = K;
478 CHKERR pipeline_mng->loopFiniteElements();
480 };
481
482 auto calculate_L_matrix = [&]() {
484
485 auto get_side_fe = [&]() {
486 auto side_fe_ptr = boost::make_shared<FaceSideEle>(mField);
487 auto det_ptr = boost::make_shared<VectorDouble>();
488 auto jac_ptr = boost::make_shared<MatrixDouble>();
489 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
490 side_fe_ptr->getOpPtrVector().push_back(
491 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
492 side_fe_ptr->getOpPtrVector().push_back(
493 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
494 side_fe_ptr->getOpPtrVector().push_back(
495 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
496 side_fe_ptr->getOpPtrVector().push_back(
497 new OpCalculateSideData("Us", "Us", FaceSideOp::OPROW));
498 side_fe_ptr->getOpPtrVector().push_back(
499 new OpCalculateSideData("Uf", "Uf", FaceSideOp::OPCOL));
500 return side_fe_ptr;
501 };
502
503 auto fe_side_ptr = boost::make_shared<SkeletonEle>(mField);
504 fe_side_ptr->getOpPtrVector().push_back(new OpLhsSkeleton(
506 matDFViscSolidPtr, boost::make_shared<Range>(rFluid),
507 boost::make_shared<Range>(rSolid), K, C, M));
508
509 auto rules = fe_side_ptr->getRuleHook = [](int, int, int p) -> int {
510 return 2 * p + 1;
511 };
512
513 CHKERR mField.loop_finite_elements(simple->getProblemName(), "INTERFACE",
514 *fe_side_ptr);
515
517 };
518
519 auto calculate_damping_matrix = [&]() {
521
522 pipeline_mng->getDomainLhsFE().reset();
523 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
524
525 CHKERR basic_ops();
526
527 pipeline_mng->getOpDomainLhsPipeline().push_back(
528 new OpK("Uf", "Uf", matDFViscFluidPtr, ptr_r_fluid));
529
530 pipeline_mng->getDomainLhsFE()->B = C;
531 CHKERR pipeline_mng->loopFiniteElements();
532
534 };
535
536 auto calculate_mass_matrix = [&]() {
538
539 pipeline_mng->getDomainLhsFE().reset();
540 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
541
542 CHKERR basic_ops();
543
544 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMass(
545 "Us", "Us",
546 [](const double, const double, const double) { return rho_solid; },
547 ptr_r_solid));
548 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMass(
549 "Uf", "Uf",
550 [](const double, const double, const double) { return rho_fluid; },
551 ptr_r_fluid));
552
553 pipeline_mng->getDomainLhsFE()->B = M;
554 CHKERR pipeline_mng->loopFiniteElements();
556 };
557
558 if (solverType == PEP_SOLVER)
559 CHKERR calculate_L_matrix();
560
561 CHKERR calculate_stiffness_matrix();
562 CHKERR calculate_damping_matrix();
563 CHKERR calculate_mass_matrix();
564
565 CHKERR MatAssemblyBegin(K, MAT_FINAL_ASSEMBLY);
566 CHKERR MatAssemblyEnd(K, MAT_FINAL_ASSEMBLY);
567 CHKERR MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
568 CHKERR MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
569 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
570 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
571
572 // auto save_range = [&](const std::string name, const Range r) {
573 // MoFEMFunctionBegin;
574 // EntityHandle out_meshset;
575 // CHKERR mField.get_moab().create_meshset(MESHSET_SET, out_meshset);
576 // CHKERR mField.get_moab().add_entities(out_meshset, r);
577 // CHKERR mField.get_moab().write_file(name.c_str(), "VTK", "",
578 // &out_meshset,
579 // 1);
580 // CHKERR mField.get_moab().delete_entities(&out_meshset, 1);
581 // MoFEMFunctionReturn(0);
582 // };
583
584 // CHKERR save_range("interface.vtk", rInterface);
585
587}
588//! [Push operators to pipeline]
589
590//! [Solve]
593 auto is_manager = mField.getInterface<ISManager>();
595
596 auto create_eps = [&](MPI_Comm comm) {
597 EPS eps;
598 CHKERR EPSCreate(comm, &eps);
599 CHKERR EPSSetOperators(eps, K, M);
600 return SmartPetscObj<EPS>(eps);
601 };
602
603 auto print_info_eps = [&](auto eps) {
605 ST st;
606 EPSType type;
607 PetscReal tol;
608 PetscInt nev, maxit, its;
609 // Optional: Get some information from the solver and display it
610 CHKERR EPSGetIterationNumber(eps, &its);
611 MOFEM_LOG_C("FSI", Sev::inform, " Number of iterations of the method: %d",
612 its);
613 CHKERR EPSGetST(eps, &st);
614 CHKERR EPSGetType(eps, &type);
615 MOFEM_LOG_C("FSI", Sev::inform, " Solution method: %s", type);
616 CHKERR EPSGetDimensions(eps, &nev, NULL, NULL);
617 MOFEM_LOG_C("FSI", Sev::inform, " Number of requested eigenvalues: %d",
618 nev);
619 CHKERR EPSGetTolerances(eps, &tol, &maxit);
620 MOFEM_LOG_C("FSI", Sev::inform, " Stopping condition: tol=%.4g, maxit=%d",
621 (double)tol, maxit);
622
623 PetscScalar eigr, eigi;
624 for (int nn = 0; nn < nev; nn++) {
625 CHKERR EPSGetEigenpair(eps, nn, &eigr, &eigi, PETSC_NULL, PETSC_NULL);
626 MOFEM_LOG_C("FSI", Sev::inform,
627 " ncov = %d eigr = %.8g eigi = %.8g (inv eigr = %.8g)", nn,
628 eigr, eigi, 1. / eigr);
629 }
630
632 };
633
634 auto setup_eps = [&](auto eps) {
636 CHKERR EPSSetFromOptions(eps);
638 };
639
640 Mat A[3] = {K.get(), C.get(), M.get()};
641
642 auto create_pep = [&](MPI_Comm comm) {
643 PEP pep;
644 CHKERR PEPCreate(comm, &pep);
645 CHKERR PEPSetOperators(pep, 3, A);
646 return SmartPetscObj<PEP>(pep);
647 };
648
649 auto setup_pep = [&](auto pep) {
651 CHKERR PEPSetFromOptions(pep);
653 };
654
655 auto print_info_pep = [&](auto pep) {
657 int nconv;
658 double kr, ki, error;
659 CHKERR PEPGetConverged(pEP, &nconv);
660 for (int j = 0; j < nconv; ++j) {
661 CHKERR PEPGetEigenpair(pep, j, &kr, &ki, PETSC_NULL, PETSC_NULL);
662 MOFEM_LOG_C("FSI", Sev::verbose, " ncov = %d eigr = %.8g eigi = %.8g", j,
663 kr, ki);
664 }
665
667 };
668
669 if (solverType == EPS_SOLVER) {
670 // Create eigensolver context
671 ePS = create_eps(mField.get_comm());
672 // Setup eps
673 CHKERR setup_eps(ePS);
674 // Solve problem
675 CHKERR EPSSolve(ePS);
676 // Print info
677 CHKERR print_info_eps(ePS);
678 } else if (solverType == PEP_SOLVER) {
679 pEP = create_pep(mField.get_comm());
680 CHKERR setup_pep(pEP);
681 CHKERR PEPSolve(pEP);
682 CHKERR print_info_pep(pEP);
683 } else {
684 SETERRQ(PETSC_COMM_SELF, MOFEM_NOT_IMPLEMENTED, "Solver not implemented");
685 }
686
688}
689//! [Solve]
690
691//! [Postprocess results]
694 auto *pipeline_mng = mField.getInterface<PipelineManager>();
695 auto *simple = mField.getInterface<Simple>();
696
697 pipeline_mng->getDomainLhsFE().reset();
698 pipeline_mng->getDomainRhsFE().reset();
699 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
700 post_proc_fe->generateReferenceElementMesh();
701
702 auto det_ptr = boost::make_shared<VectorDouble>();
703 auto jac_ptr = boost::make_shared<MatrixDouble>();
704 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
705 post_proc_fe->getOpPtrVector().push_back(
706 new OpCalculateHOJac<SPACE_DIM>(jac_ptr));
707 post_proc_fe->getOpPtrVector().push_back(
708 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
709 post_proc_fe->getOpPtrVector().push_back(
710 new OpSetHOInvJacToScalarBases<SPACE_DIM>(H1, inv_jac_ptr));
711
712 auto us_ptr = boost::make_shared<MatrixDouble>();
713 auto grad_ptr = boost::make_shared<MatrixDouble>();
714 auto strain_s_ptr = boost::make_shared<MatrixDouble>();
715 auto uf_ptr = boost::make_shared<MatrixDouble>();
716 auto strain_f_ptr = boost::make_shared<MatrixDouble>();
717
718 // auto stress_ptr = boost::make_shared<MatrixDouble>();
719
720 post_proc_fe->getOpPtrVector().push_back(
722 post_proc_fe->getOpPtrVector().push_back(
724 post_proc_fe->getOpPtrVector().push_back(
725 new OpSymmetrizeTensor<SPACE_DIM>("Us", grad_ptr, strain_s_ptr));
726
727 post_proc_fe->getOpPtrVector().push_back(
729 post_proc_fe->getOpPtrVector().push_back(
731 post_proc_fe->getOpPtrVector().push_back(
732 new OpSymmetrizeTensor<SPACE_DIM>("Uf", grad_ptr, strain_f_ptr));
733
735
736 post_proc_fe->getOpPtrVector().push_back(
737
738 new OpPPMap(post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts,
739
741
742 OpPPMap::DataMapMat{{"Us", us_ptr}, {"Uf", uf_ptr}},
743
745
746 OpPPMap::DataMapMat{{"STRAINs", strain_s_ptr},
747 {"STRAINf", strain_f_ptr}}
748
749 )
750
751 );
752
753 auto bc_mng = mField.getInterface<BcManager>();
754 auto &bc_map = bc_mng->getBcMapByBlockName();
755
756 auto op_tag_handle = new DomainEleOp(NOSPACE, DomainEleOp::OPSPACE);
757 op_tag_handle->doWorkRhsHook = [&](DataOperator *op_ptr, int side,
761
762 auto get_tag = [&]() {
763 int def_val = -1;
764 Tag th;
765 CHK_MOAB_THROW(post_proc_fe->postProcMesh.tag_get_handle(
766 "BLOCK", 1, MB_TYPE_INTEGER, th,
767 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val),
768 "Can not create tag on post-proc mesh");
769 return th;
770 };
771
772 auto fe_ent = static_cast<DomainEleOp *>(op_ptr)->getFEEntityHandle();
773
774 int marker = 0;
775 if (rFluid.find(fe_ent) != rFluid.end()) {
776 marker += 1;
777 }
778 if (rSolid.find(fe_ent) != rSolid.end()) {
779 marker += 2;
780 }
781
782 auto &nodes_vec = post_proc_fe->mapGaussPts;
784 post_proc_fe->postProcMesh.tag_clear_data(
785 get_tag(), &*nodes_vec.begin(), nodes_vec.size(), &marker),
786 "Can not set tag");
787
789 };
790
791 post_proc_fe->getOpPtrVector().push_back(op_tag_handle);
792
793 pipeline_mng->getDomainRhsFE() = post_proc_fe;
794
795 auto dm = simple->getDM();
796
797 auto save_vec = [&](auto file_name, auto v) {
799 CHKERR DMoFEMMeshToLocalVector(dm, v, INSERT_VALUES, SCATTER_REVERSE);
800 CHKERR pipeline_mng->loopFiniteElements();
801 post_proc_fe->writeFile(file_name);
803 };
804
805 auto post_proc_eps = [&](auto eps) {
807 auto v = smartCreateDMVector(dm);
808
809 PetscInt nev;
810 CHKERR EPSGetDimensions(ePS, &nev, NULL, NULL);
811 PetscScalar eigr, eigi, nrm2r;
812 for (int nn = 0; nn < nev; nn++) {
813 CHKERR EPSGetEigenpair(ePS, nn, &eigr, &eigi, v, PETSC_NULL);
814 CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
815 CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
816 CHKERR VecNorm(v, NORM_2, &nrm2r);
818 "FSI", Sev::verbose,
819 " ncov = %d omega2 = %.8g omega = %.8g osc frequency = %.8g kHz",
820 nn, eigr, std::sqrt(std::abs(eigr)),
821 std::sqrt(std::abs(eigr)) / (2 * M_PI));
822
823 CHKERR save_vec(
824 "out_eig_" + boost::lexical_cast<std::string>(nn) + ".h5m", v);
825 }
826
828 };
829
830 auto post_proc_pep = [&](auto eps) {
832
833 auto re_vec = smartCreateDMVector(dm);
834 auto im_vec = smartVectorDuplicate(re_vec);
835
836 PetscInt nev;
837 CHKERR PEPGetConverged(pEP, &nev);
838 PetscScalar ki, kr, error, nrm2r;
839
840 int idx = 0; // counting from zero, increment first before log
841 for (int nn = 0; nn < nev; nn++) {
842
843 CHKERR PEPGetEigenpair(pEP, nn, &kr, &ki, re_vec, im_vec);
844 CHKERR PEPComputeError(pEP, nn, PEP_ERROR_BACKWARD, &error);
845
846 CHKERR VecGhostUpdateBegin(re_vec, INSERT_VALUES, SCATTER_FORWARD);
847 CHKERR VecGhostUpdateEnd(re_vec, INSERT_VALUES, SCATTER_FORWARD);
848 CHKERR VecGhostUpdateBegin(im_vec, INSERT_VALUES, SCATTER_FORWARD);
849 CHKERR VecGhostUpdateEnd(im_vec, INSERT_VALUES, SCATTER_FORWARD);
850
851 if (std::abs(ki) > eig_threshold &&
852 kr < std::numeric_limits<float>::epsilon()) {
853 CHKERR save_vec("out_pep_re_" + boost::lexical_cast<std::string>(idx) +
854 ".h5m",
855 re_vec);
856 CHKERR save_vec("out_pep_im_" + boost::lexical_cast<std::string>(idx) +
857 ".h5m",
858 im_vec);
859 ++idx;
860 }
861
862 MOFEM_LOG_C("FSI", Sev::inform,
863 "idx = %d ncov = %d eigr = %.8g eigi = %.8g angular "
864 "frequency = %.8g kHz "
865 "error %.8g",
866 idx - 1, nn, kr, ki, ki / (2 * M_PI), error);
867 }
869 };
870
871 if (solverType == EPS_SOLVER)
872 CHKERR post_proc_eps(ePS);
873 else if (solverType == PEP_SOLVER)
874 CHKERR post_proc_pep(pEP);
875
877}
878//! [Postprocess results]
879
880static char help[] = "...\n\n";
881
882int main(int argc, char *argv[]) {
883
884 // Initialisation of MoFEM/PETSc and MOAB data structures
885 const char param_file[] = "param_file.petsc";
886 SlepcInitialize(&argc, &argv, param_file, help);
888
889 // Add logging channel for FSI
890 auto core_log = logging::core::get();
891 core_log->add_sink(LogManager::createSink(LogManager::getStrmWorld(), "FSI"));
892 LogManager::setLog("FSI");
893 MOFEM_LOG_TAG("FSI", "FSI");
894
895 try {
896
897 //! [Register MoFEM discrete manager in PETSc]
898 DMType dm_name = "DMMOFEM";
899 CHKERR DMRegister_MoFEM(dm_name);
900 //! [Register MoFEM discrete manager in PETSc
901
902 //! [Create MoAB]
903 moab::Core mb_instance; ///< mesh database
904 moab::Interface &moab = mb_instance; ///< mesh database interface
905 //! [Create MoAB]
906
907 //! [Create MoFEM]
908 MoFEM::Core core(moab); ///< finite element database
909 MoFEM::Interface &m_field = core; ///< finite element database insterface
910 //! [Create MoFEM]
911
912 //! [FSI]
913 FSI ex(m_field);
914 CHKERR ex.runProblem();
915 //! [FSI]
916 }
918
919 SlepcFinalize();
921}
922
923// data for skeleton computation
924std::array<std::vector<VectorInt>, 2> indicesUsSideMap;
925std::array<std::vector<VectorInt>, 2> indicesUfSideMap;
926
927std::array<std::vector<MatrixDouble>, 2> baseUfSideMap;
928std::array<std::vector<MatrixDouble>, 2> diffUfBaseSideMap;
929std::array<std::vector<MatrixDouble>, 2> baseUsSideMap;
930std::array<std::vector<MatrixDouble>, 2> diffUsBaseSideMap;
931
932std::array<EntityHandle, 2> sideHandle;
933std::array<double, 2> areaMap; // area/volume of elements on the side
934std::array<int, 2> senseMap; // orientation of local element edge/face in
935 // respect to global orientation of edge/face
936
938 std::string col_field_name,
939 UserDataOperator::OpType type)
940 : FaceSideOp(row_field_name, col_field_name, type) {}
941
943 EntData &data) {
945
946 const auto nb_in_loop = getFEMethod()->nInTheLoop;
947
948 auto clear = [](auto nb) {
949 indicesUsSideMap[nb].clear();
950 indicesUfSideMap[nb].clear();
951 baseUsSideMap[nb].clear();
952 diffUsBaseSideMap[nb].clear();
953 baseUfSideMap[nb].clear();
954 diffUfBaseSideMap[nb].clear();
955 };
956
957 if (getOpType() == OPROW) {
958
959 if (type == MBVERTEX) {
960 areaMap[nb_in_loop] = getMeasure();
961 senseMap[nb_in_loop] = getSkeletonSense();
962 sideHandle[nb_in_loop] = getFEEntityHandle();
963 if (!nb_in_loop) {
964 clear(0);
965 clear(1);
966 sideHandle[1] = 0;
967 areaMap[1] = 0;
968 senseMap[1] = 0;
969 }
970 }
971
972 const auto nb_dofs = data.getIndices().size();
973 if (nb_dofs) {
974 indicesUsSideMap[nb_in_loop].push_back(data.getIndices());
975 baseUsSideMap[nb_in_loop].push_back(
976 data.getN(BaseDerivatives::ZeroDerivative));
977 diffUsBaseSideMap[nb_in_loop].push_back(
978 data.getN(BaseDerivatives::FirstDerivative));
979 }
980 } else if (getOpType() == OPCOL) {
981
982 if (type == MBVERTEX) {
983 if (areaMap[nb_in_loop] != getMeasure()) {
984 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
985 "Should be the same");
986 }
987 if (senseMap[nb_in_loop] != getSkeletonSense()) {
988 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
989 "Should be the same");
990 }
991 if (sideHandle[nb_in_loop] != getFEEntityHandle()) {
992 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
993 "Should be the same");
994 };
995 }
996
997 const auto nb_dofs = data.getIndices().size();
998 if (nb_dofs) {
999 indicesUfSideMap[nb_in_loop].push_back(data.getIndices());
1000 baseUfSideMap[nb_in_loop].push_back(
1001 data.getN(BaseDerivatives::ZeroDerivative));
1002 diffUfBaseSideMap[nb_in_loop].push_back(
1003 data.getN(BaseDerivatives::FirstDerivative));
1004 }
1005
1006 } else {
1007 SETERRQ1(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Wrong OpType %s",
1008 OpTypeNames[getOpType()]);
1009 }
1010
1012}
1013
1015
1016template <typename T> inline auto get_ntensor(T &base_mat) {
1018 &*base_mat.data().begin());
1019};
1020
1021template <typename T> inline auto get_ntensor(T &base_mat, int gg, int bb) {
1022 double *ptr = &base_mat(gg, bb);
1024};
1025
1026template <typename T> inline auto get_diff_ntensor(T &base_mat) {
1027 double *ptr = &*base_mat.data().begin();
1028 return getFTensor1FromPtr<SPACE_DIM>(ptr);
1029};
1030
1031template <typename T>
1032inline auto get_diff_ntensor(T &base_mat, int gg, int bb) {
1033 double *ptr = &base_mat(gg, SPACE_DIM * bb);
1034 return getFTensor1FromPtr<SPACE_DIM>(ptr);
1035};
1036
1037/**
1038 * @brief Construct a new OpH1LhsSkeleton
1039 *
1040 * @param side_fe_ptr pointer to FE to evaluate side elements
1041 */
1043 boost::shared_ptr<FaceSideEle> side_fe_ptr,
1044 boost::shared_ptr<MatrixDouble> mat_d_fluid_ptr,
1045 boost::shared_ptr<MatrixDouble> mat_d_solid_ptr,
1046 boost::shared_ptr<MatrixDouble> mat_d_visc_fluid_ptr,
1047 boost::shared_ptr<MatrixDouble> mat_d_visc_solid_ptr,
1048 boost::shared_ptr<Range> r_fluid_ptr, boost::shared_ptr<Range> r_solid_ptr,
1050 : SkeletonEleOp(NOSPACE, SkeletonEleOp::OPSPACE), sideFEPtr(side_fe_ptr),
1051 dMatFluidPtr(mat_d_fluid_ptr), dMatSolidPtr(mat_d_solid_ptr),
1052 dViscFluidPtr(mat_d_visc_fluid_ptr), dViscSolidPtr(mat_d_visc_solid_ptr),
1053 rFluidPtr(r_fluid_ptr), rSolidPtr(r_solid_ptr), K(K), C(C), M(M) {}
1054
1058
1059 // get normal of the face or edge
1060 auto t_normal = getFTensor1Normal();
1061 t_normal(i) /= sqrt(t_normal(j) * t_normal(j));
1062
1063 // Collect data from side domian elements
1064 CHKERR loopSide("dFE", sideFEPtr.get(), SPACE_DIM);
1065 const auto in_the_loop =
1066 sideFEPtr->nInTheLoop; // return number of elements on the side
1067
1068 std::array<boost::shared_ptr<MatrixDouble>, 2> d_ptr;
1069 std::array<boost::shared_ptr<MatrixDouble>, 2> d_vis_ptr;
1070
1071 std::array<std::vector<VectorInt>, 2> indices_side_map;
1072 std::array<std::vector<MatrixDouble>, 2> base_side_map;
1073 std::array<std::vector<MatrixDouble>, 2> diff_base_side_map;
1074
1075 // MOFEM_LOG("SELF", Sev::inform)
1076 // << "AAA " << sideHandle[0] << " " << sideHandle[1];
1077
1078 // auto save_range = [&](const std::string name, const Range r) {
1079 // MoFEMFunctionBegin;
1080 // EntityHandle out_meshset;
1081 // CHKERR getPtrFE()->mField.get_moab().create_meshset(MESHSET_SET,
1082 // out_meshset);
1083 // CHKERR getPtrFE()->mField.get_moab().add_entities(out_meshset, r);
1084 // CHKERR getPtrFE()->mField.get_moab().write_file(name.c_str(), "VTK", "",
1085 // &out_meshset, 1);
1086 // CHKERR getPtrFE()->mField.get_moab().delete_entities(&out_meshset, 1);
1087 // MoFEMFunctionReturn(0);
1088 // };
1089
1090 // Tag th_n;
1091 // double def[] = {0, 0, 0};
1092 // CHKERR getPtrFE()->mField.get_moab().tag_get_handle(
1093 // "NORMAL", 3, MB_TYPE_DOUBLE, th_n, MB_TAG_CREAT | MB_TAG_SPARSE, def);
1094
1095 auto set_row_col = [&](auto a, auto b) {
1096 d_ptr[a] = dMatSolidPtr;
1097 d_ptr[b] = dMatFluidPtr;
1098 d_vis_ptr[a] = dViscSolidPtr;
1099 d_vis_ptr[b] = dViscFluidPtr;
1100 indices_side_map[a].swap(indicesUsSideMap[a]);
1101 indices_side_map[b].swap(indicesUfSideMap[b]);
1102 base_side_map[a].swap(baseUsSideMap[a]);
1103 base_side_map[b].swap(baseUfSideMap[b]);
1104 diff_base_side_map[a].swap(diffUsBaseSideMap[a]);
1105 diff_base_side_map[b].swap(diffUfBaseSideMap[b]);
1106 };
1107
1108 if (rSolidPtr->find(sideHandle[0]) != rSolidPtr->end()) {
1109
1110 if (rFluidPtr->find(sideHandle[1]) == rFluidPtr->end()) {
1111 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1112 "Other side should be fluid");
1113 }
1114
1115 set_row_col(0, 1);
1116
1117 if (senseMap[0] < 0)
1118 t_normal(i) *= -1;
1119
1120 } else if (rFluidPtr->find(sideHandle[0]) != rFluidPtr->end()) {
1121
1122 if (rSolidPtr->find(sideHandle[1]) == rSolidPtr->end()) {
1123 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1124 "Other side should be solid");
1125 }
1126
1127 set_row_col(1, 0);
1128
1129 if (senseMap[1] < 0)
1130 t_normal(i) *= -1;
1131
1132 } else {
1133 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1134 "Side is not fluid or solid");
1135 }
1136
1137 // auto fe_h = getFEEntityHandle();
1138 // CHKERR getPtrFE()->mField.get_moab().tag_set_data(th_n, &fe_h, 1,
1139 // &t_normal(0));
1140
1141 // calculate penalty
1142 const double s = getMeasure() / (areaMap[0] + areaMap[1]);
1143 const double p = pow(penalty, order) * s;
1144
1145 // get number of integration points on face
1146 const size_t nb_integration_pts = getGaussPts().size2();
1147
1148 // when boundary edge/face is evaluated, and 2 when is skeleton edge/face.
1149 const double beta = 1. / (in_the_loop + 1);
1150
1151 auto integrate = [&](auto sense_row, auto &row_ind, auto &row, auto &row_diff,
1152 auto sense_col, auto &col_ind, auto &col, auto &col_diff,
1153 auto &d_row, auto &d_col, auto &d_vis_row,
1154 auto &d_vis_col) {
1156
1157 // number of dofs, for homogenous approximation this value is
1158 // constant.
1159 const auto nb_rows = row_ind.size();
1160 const auto nb_cols = col_ind.size();
1161
1162 const auto nb_row_base_functions = row.size2();
1163
1164 auto t_D_row = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_row);
1165 auto t_D_col = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_col);
1166
1167 auto t_vis_D_row =
1168 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_vis_row);
1169 auto t_vis_D_col =
1170 getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(d_vis_col);
1171
1172 if (nb_cols) {
1173
1174 // resize local element matrix
1175 locMatK.resize(nb_rows, nb_cols, false);
1176 locMatK.clear();
1177 locMatC.resize(nb_rows, nb_cols, false);
1178 locMatC.clear();
1179 locMatM.resize(nb_rows, nb_cols, false);
1180 locMatM.clear();
1181
1182 // get base functions, and integration weights
1183 auto t_row_base = get_ntensor(row);
1184 auto t_diff_row_base = get_diff_ntensor(row_diff);
1185
1186 constexpr auto t_kd = FTensor::Kronecker_Delta<double>();
1188 t_P(i, j) = t_normal(i) * t_normal(j);
1189
1190 auto t_w = getFTensor0IntegrationWeight();
1191
1192 // iterate integration points on face/edge
1193 for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
1194
1195 // t_w is integration weight, and measure is element area, or
1196 // volume, depending if problem is in 2d/3d.
1197 const double alpha = getMeasure() * t_w;
1198
1199 // iterate rows
1200 size_t rr = 0;
1201 for (; rr != nb_rows / SPACE_DIM; ++rr) {
1202
1203 auto get_vn_plus = [&](auto &t_D) {
1205 t_vn_plus(i, j, k) =
1206 (beta * phi / p) *
1207 (t_kd(i, m) * t_D(m, j, k, l) * t_diff_row_base(l));
1208 return t_vn_plus;
1209 };
1210
1211 auto get_vn = [&]() {
1213 t_vn(i, j, k) = t_normal(j) * (t_kd(i, k) * t_row_base * sense_row);
1214 return t_vn;
1215 };
1216
1217 auto t_vn_plus = get_vn_plus(t_D_row);
1218 auto t_vn = get_vn();
1219 auto t_vn_lambda_plus = get_vn_plus(t_vis_D_row);
1220
1221 // get base functions on columns
1222 auto t_col_base = get_ntensor(col, gg, 0);
1223 auto t_diff_col_base = get_diff_ntensor(col_diff, gg, 0);
1224
1225 auto t_mat_K = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1226 locMatK, SPACE_DIM * rr);
1227 auto t_mat_C = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1228 locMatC, SPACE_DIM * rr);
1229 auto t_mat_M = getFTensor2FromArray<SPACE_DIM, SPACE_DIM, SPACE_DIM>(
1230 locMatM, SPACE_DIM * rr);
1231
1232 // iterate columns
1233 for (size_t cc = 0; cc != nb_cols / SPACE_DIM; ++cc) {
1234
1235 auto get_un_plus = [&](auto &t_D) {
1237 t_un_plus;
1238 t_un_plus(i, j, k) =
1239 beta * ((t_P(i, m) * t_D(m, j, k, l)) * t_diff_col_base(l));
1240 return t_un_plus;
1241 };
1242
1243 auto get_un = [&]() {
1245 t_un(i, j, k) =
1246 (t_normal(j) * (t_P(i, k) * t_col_base * sense_col));
1247 return t_un;
1248 };
1249
1250 auto t_un_plus = get_un_plus(t_D_col);
1251 auto t_un = get_un();
1252 auto t_un_lambda_plus = get_un_plus(t_vis_D_col);
1253
1254 // K
1255 t_mat_K(i, j) -=
1256 alpha * (t_vn(m, n, i) - t_vn_plus(m, n, i)) *
1257 ((-p) * (t_un(m, n, j) - (t_un_plus(m, n, j) / p)));
1258 t_mat_K(i, j) -= alpha * (t_vn_plus(m, n, i) * t_un_plus(m, n, j));
1259
1260 // C1
1261 t_mat_C(i, j) -=
1262 alpha * (-t_vn_lambda_plus(m, n, i)) *
1263 ((-p) * (t_un(m, n, j) - (t_un_plus(m, n, j) / p)));
1264 t_mat_C(i, j) -=
1265 alpha * (t_vn_lambda_plus(m, n, i) * t_un_plus(m, n, j));
1266
1267 // C2
1268 t_mat_C(i, j) -= alpha * (t_vn(m, n, i) - t_vn_plus(m, n, i)) *
1269 ((-p) * (t_un_lambda_plus(m, n, j) / (-p)));
1270 t_mat_C(i, j) -=
1271 alpha * (t_vn_plus(m, n, i) * t_un_lambda_plus(m, n, j));
1272
1273 // M
1274 t_mat_M(i, j) -= alpha * (-t_vn_lambda_plus(m, n, i)) *
1275 ((-p) * (t_un_lambda_plus(m, n, j) / (-p)));
1276 t_mat_M(i, j) -=
1277 alpha * (t_vn_lambda_plus(m, n, i) * t_un_lambda_plus(m, n, j));
1278
1279 // move to next column base and element of matrix
1280 ++t_col_base;
1281 ++t_diff_col_base;
1282 ++t_mat_K;
1283 ++t_mat_C;
1284 ++t_mat_M;
1285 }
1286
1287 // move to next row base
1288 ++t_row_base;
1289 ++t_diff_row_base;
1290 }
1291
1292 for (; rr < nb_row_base_functions; ++rr) {
1293 ++t_row_base;
1294 ++t_diff_row_base;
1295 }
1296
1297 ++t_w;
1298 }
1299
1300 // assemble system
1301
1302 CHKERR ::MatSetValues(K, nb_rows, &*row_ind.begin(), col_ind.size(),
1303 &*col_ind.begin(), &*locMatK.data().begin(),
1304 ADD_VALUES);
1305 CHKERR ::MatSetValues(C, nb_rows, &*row_ind.begin(), col_ind.size(),
1306 &*col_ind.begin(), &*locMatC.data().begin(),
1307 ADD_VALUES);
1308 CHKERR ::MatSetValues(M, nb_rows, &*row_ind.begin(), col_ind.size(),
1309 &*col_ind.begin(), &*locMatM.data().begin(),
1310 ADD_VALUES);
1311 }
1312
1314 };
1315
1316 // iterate the sides rows
1317 for (auto s0 : {LEFT_SIDE, RIGHT_SIDE}) {
1318 const auto sense_row = senseMap[s0];
1319
1320 for (auto x0 = 0; x0 != indices_side_map[s0].size(); ++x0) {
1321
1322 for (auto s1 : {LEFT_SIDE, RIGHT_SIDE}) {
1323 const auto sense_col = senseMap[s1];
1324
1325 for (auto x1 = 0; x1 != indices_side_map[s1].size(); ++x1) {
1326
1327 CHKERR integrate(sense_row, indices_side_map[s0][x0],
1328 base_side_map[s0][x0], diff_base_side_map[s0][x0],
1329
1330 sense_col, indices_side_map[s1][x1],
1331 base_side_map[s1][x1], diff_base_side_map[s1][x1],
1332
1333 *d_ptr[s0], *d_ptr[s1], *d_vis_ptr[s0],
1334 *d_vis_ptr[s1]
1335
1336 );
1337 }
1338 }
1339 }
1340 }
1341
1343}
1344
static Index< 'M', 3 > M
static Index< 'p', 3 > p
std::string param_file
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
int main()
constexpr double a
static const double eps
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
Kronecker Delta class.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
@ H1
continuous field
Definition definitions.h:85
@ NOSPACE
Definition definitions.h:83
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
@ BLOCKSET
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ MOFEM_NOT_IMPLEMENTED
Definition definitions.h:32
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
const int dim
#define MU(E, NU)
Definition fem_tools.h:23
auto get_diff_ntensor(T &base_mat)
std::array< std::vector< MatrixDouble >, 2 > baseUfSideMap
ElementsAndOps< SPACE_DIM >::DomainEleOp DomainEleOp
FTensor::Index< 'j', SPACE_DIM > j
static char help[]
[Postprocess results]
FTensor::Index< 'k', SPACE_DIM > k
std::array< std::vector< MatrixDouble >, 2 > baseUsSideMap
FTensor::Index< 'i', SPACE_DIM > i
auto get_ntensor(T &base_mat)
ElementsAndOps< SPACE_DIM >::SkeletonEle SkeletonEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
constexpr int SPACE_DIM
FTensor::Index< 'l', SPACE_DIM > l
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
std::array< double, 2 > areaMap
std::array< int, 2 > senseMap
std::array< std::vector< VectorInt >, 2 > indicesUfSideMap
FTensor::Index< 'n', SPACE_DIM > n
std::array< EntityHandle, 2 > sideHandle
std::array< std::vector< MatrixDouble >, 2 > diffUfBaseSideMap
static double penalty
FTensor::Index< 'm', SPACE_DIM > m
ElementsAndOps< SPACE_DIM >::FaceSideEle FaceSideEle
std::array< std::vector< MatrixDouble >, 2 > diffUsBaseSideMap
static double phi
std::array< std::vector< VectorInt >, 2 > indicesUsSideMap
auto integration_rule
constexpr auto t_kd
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:509
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition DMMoFEM.cpp:1183
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:47
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_filed)=0
set finite element field data
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode loop_finite_elements(const std::string problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
auto marker
set bit to marker
const double v
phase velocity of light in medium (cm/ns)
constexpr IntegrationType G
Definition level_set.cpp:33
double tol
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
DEPRECATED auto smartCreateDMVector(DM dm)
Definition DMMoFEM.hpp:1011
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DEPRECATED SmartPetscObj< Mat > smartMatDuplicate(Mat mat, MatDuplicateOption op)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition plastic.cpp:13
constexpr auto size_symm
Definition plastic.cpp:42
std::array< double, 2 > areaMap
Definition plate.cpp:102
std::array< int, 2 > senseMap
Definition plate.cpp:103
constexpr auto field_name
static constexpr int approx_order
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
[Read mesh]
SmartPetscObj< Mat > M
SmartPetscObj< Mat > C
MoFEMErrorCode assembleSystem()
[Boundary condition]
boost::shared_ptr< MatrixDouble > matDFluidPtr
MoFEMErrorCode createCommonData()
[Create common data]
boost::shared_ptr< MatrixDouble > matDFViscSolidPtr
MoFEMErrorCode runProblem()
[Create common data]
SmartPetscObj< EPS > ePS
boost::shared_ptr< MatrixDouble > matDFViscFluidPtr
boost::shared_ptr< MatrixDouble > matStressPtr
boost::shared_ptr< MatrixDouble > matDSolidPtr
boost::shared_ptr< MatrixDouble > matGradPtr
SmartPetscObj< PEP > pEP
FSI(MoFEM::Interface &m_field)
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEMErrorCode boundaryCondition()
[Set up problem]
SmartPetscObj< Mat > K
boost::shared_ptr< MatrixDouble > matStrainPtr
MoFEMErrorCode outputResults()
[Solve]
MoFEMErrorCode readMesh()
[Run problem]
Simple interface for fast problem set-up.
Definition BcManager.hpp:25
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:112
base operator to do operations at Gauss Pt. level
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
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.
Base face element used to integrate on skeleton.
@ OPSPACE
operator do Work is execute on space data
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
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.
Post post-proc data at points from hash maps.
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator tp collect data from elements on the side of Edge/Face.
OpCalculateSideData(std::string field_name, std::string col_field_name, UserDataOperator::OpType type)
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator the left hand side matrix.
boost::shared_ptr< MatrixDouble > dMatSolidPtr
boost::shared_ptr< Range > rSolidPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpLhsSkeleton(boost::shared_ptr< FaceSideEle > side_fe_ptr, boost::shared_ptr< MatrixDouble > d_mat_fluid_ptr, boost::shared_ptr< MatrixDouble > d_mat_solid_ptr, boost::shared_ptr< MatrixDouble > d_mat_visc_fluid_ptr, boost::shared_ptr< MatrixDouble > d_mat_visc_solid_ptr, boost::shared_ptr< Range > r_visc_fluid_ptr, boost::shared_ptr< Range > r_visc_solid_ptr, SmartPetscObj< Mat > K, SmartPetscObj< Mat > C, SmartPetscObj< Mat > M)
Construct a new OpH1LhsSkeleton.
boost::shared_ptr< MatrixDouble > dMatFluidPtr
boost::shared_ptr< Range > rFluidPtr
boost::shared_ptr< MatrixDouble > dViscSolidPtr
boost::shared_ptr< FaceSideEle > sideFEPtr
pointer to element to get data on edge/face sides
boost::shared_ptr< MatrixDouble > dViscFluidPtr
Post post-proc data at points from hash maps.
Postprocess on face.
Post processing.