v0.13.1
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]
282}
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,...)
Definition: LogManager.hpp:304
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
static const double eps
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
Kronecker Delta class.
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ 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()
Definition: definitions.h:447
@ 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 ...
Definition: definitions.h:346
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
@ BLOCKSET
Definition: definitions.h:148
@ 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()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
double eig_threshold
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
constexpr int SPACE_DIM
FTensor::Index< 'l', SPACE_DIM > l
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
std::array< double, 2 > areaMap
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
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: DMMMoFEM.cpp:470
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMMoFEM.cpp:1144
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
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
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_data(const std::string &fe_name, const std::string &name_filed)=0
set finite element field data
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
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
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)
double tol
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SmartPetscObj< Mat > smartMatDuplicate(Mat &mat, MatDuplicateOption op)
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:10
constexpr auto size_symm
Definition: plastic.cpp:33
constexpr AssemblyType A
Definition: plastic.cpp:35
constexpr IntegrationType G
Definition: plastic.cpp:36
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:23
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:165
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.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
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:26
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Base volume element used to integrate on skeleton.
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.