v0.14.0
Loading...
Searching...
No Matches
Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
FSI Struct Reference
Collaboration diagram for FSI:
[legend]

Public Member Functions

 FSI (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 [Create common data] More...
 

Private Types

enum  SolverType { EPS_SOLVER , PEP_SOLVER , LASTSOLVER }
 

Private Member Functions

MoFEMErrorCode readMesh ()
 [Run problem] More...
 
MoFEMErrorCode setupProblem ()
 [Read mesh] More...
 
MoFEMErrorCode createCommonData ()
 [Create common data] More...
 
MoFEMErrorCode boundaryCondition ()
 [Set up problem] More...
 
MoFEMErrorCode assembleSystem ()
 [Boundary condition] More...
 
MoFEMErrorCode solveSystem ()
 [Push operators to pipeline] More...
 
MoFEMErrorCode outputResults ()
 [Solve] More...
 

Private Attributes

MoFEM::InterfacemField
 
boost::shared_ptr< MatrixDouble > matGradPtr
 
boost::shared_ptr< MatrixDouble > matStrainPtr
 
boost::shared_ptr< MatrixDouble > matStressPtr
 
boost::shared_ptr< MatrixDouble > matDSolidPtr
 
boost::shared_ptr< MatrixDouble > matDFluidPtr
 
boost::shared_ptr< MatrixDouble > matDFViscFluidPtr
 
boost::shared_ptr< MatrixDouble > matDFViscSolidPtr
 
SmartPetscObj< Mat > M
 
SmartPetscObj< Mat > C
 
SmartPetscObj< Mat > K
 
SmartPetscObj< EPS > ePS
 
SmartPetscObj< PEP > pEP
 
int solverType
 
Range rSolid
 
Range rFluid
 
Range rInterface
 

Detailed Description

Definition at line 146 of file fluid_structure_eigenproblem.cpp.

Member Enumeration Documentation

◆ SolverType

enum FSI::SolverType
private
Enumerator
EPS_SOLVER 
PEP_SOLVER 
LASTSOLVER 

Definition at line 178 of file fluid_structure_eigenproblem.cpp.

Constructor & Destructor Documentation

◆ FSI()

FSI::FSI ( MoFEM::Interface m_field)
inline

Definition at line 148 of file fluid_structure_eigenproblem.cpp.

148: mField(m_field) {}
MoFEM::Interface & mField

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode FSI::assembleSystem ( )
private

[Boundary condition]

[Push operators to pipeline]

Definition at line 424 of file fluid_structure_eigenproblem.cpp.

424 {
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}
static Index< 'p', 3 > p
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr int SPACE_DIM
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
auto integration_rule
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition: DMMoFEM.cpp:1183
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.
DEPRECATED SmartPetscObj< Mat > smartMatDuplicate(Mat mat, MatDuplicateOption op)
static constexpr int approx_order
SmartPetscObj< Mat > M
SmartPetscObj< Mat > C
boost::shared_ptr< MatrixDouble > matDFluidPtr
boost::shared_ptr< MatrixDouble > matDFViscSolidPtr
boost::shared_ptr< MatrixDouble > matDFViscFluidPtr
boost::shared_ptr< MatrixDouble > matDSolidPtr
SmartPetscObj< Mat > K
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Operator tp collect data from elements on the side of Edge/Face.
Operator the left hand side matrix.

◆ boundaryCondition()

MoFEMErrorCode FSI::boundaryCondition ( )
private

[Set up problem]

[Boundary condition]

Definition at line 393 of file fluid_structure_eigenproblem.cpp.

393 {
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}

◆ createCommonData()

MoFEMErrorCode FSI::createCommonData ( )
private

[Create common data]

Definition at line 185 of file fluid_structure_eigenproblem.cpp.

185 {
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}
Kronecker Delta class symmetric.
#define MU(E, NU)
Definition: fem_tools.h:23
FTensor::Index< 'j', SPACE_DIM > j
FTensor::Index< 'k', SPACE_DIM > k
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', SPACE_DIM > l
static double penalty
static double phi
constexpr auto t_kd
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
constexpr AssemblyType A
constexpr auto size_symm
Definition: plastic.cpp:42
boost::shared_ptr< MatrixDouble > matStressPtr
boost::shared_ptr< MatrixDouble > matGradPtr
boost::shared_ptr< MatrixDouble > matStrainPtr

◆ outputResults()

MoFEMErrorCode FSI::outputResults ( )
private

[Solve]

[Postprocess results]

Definition at line 692 of file fluid_structure_eigenproblem.cpp.

692 {
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}
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
static const double eps
@ NOSPACE
Definition: definitions.h:83
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:576
ElementsAndOps< SPACE_DIM >::DomainEleOp DomainEleOp
double eig_threshold
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
auto marker
set bit to marker
const double v
phase velocity of light in medium (cm/ns)
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
DEPRECATED auto smartCreateDMVector(DM dm)
Definition: DMMoFEM.hpp:1011
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
SmartPetscObj< EPS > ePS
SmartPetscObj< PEP > pEP
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
Definition: BcManager.hpp:200
base operator to do operations at Gauss Pt. level
Data on single entity (This is passed as argument to DataOperator::doWork)
@ OPSPACE
operator do Work is execute on space data
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
Post post-proc data at points from hash maps.

◆ readMesh()

MoFEMErrorCode FSI::readMesh ( )
private

[Run problem]

[Read mesh]

Definition at line 286 of file fluid_structure_eigenproblem.cpp.

286 {
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}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ BLOCKSET
Definition: definitions.h:148
#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 _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.
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
int r
Definition: sdf.py:8
#define EXECUTABLE_DIMENSION
Definition: plastic.cpp:13
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0

◆ runProblem()

MoFEMErrorCode FSI::runProblem ( )

[Create common data]

[Run problem]

Definition at line 272 of file fluid_structure_eigenproblem.cpp.

272 {
282}
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode outputResults()
[Solve]
MoFEMErrorCode readMesh()
[Run problem]

◆ setupProblem()

MoFEMErrorCode FSI::setupProblem ( )
private

[Read mesh]

[Set up problem]

Definition at line 348 of file fluid_structure_eigenproblem.cpp.

348 {
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}
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
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
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Problem manager is used to build and partition problems.

◆ solveSystem()

MoFEMErrorCode FSI::solveSystem ( )
private

[Push operators to pipeline]

[Solve]

Definition at line 591 of file fluid_structure_eigenproblem.cpp.

591 {
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}
@ MOFEM_NOT_IMPLEMENTED
Definition: definitions.h:32
double tol
virtual MPI_Comm & get_comm() const =0
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
intrusive_ptr for managing petsc objects

Member Data Documentation

◆ C

SmartPetscObj<Mat> FSI::C
private

Definition at line 172 of file fluid_structure_eigenproblem.cpp.

◆ ePS

SmartPetscObj<EPS> FSI::ePS
private

Definition at line 175 of file fluid_structure_eigenproblem.cpp.

◆ K

SmartPetscObj<Mat> FSI::K
private

Definition at line 173 of file fluid_structure_eigenproblem.cpp.

◆ M

SmartPetscObj<Mat> FSI::M
private

Definition at line 171 of file fluid_structure_eigenproblem.cpp.

◆ matDFluidPtr

boost::shared_ptr<MatrixDouble> FSI::matDFluidPtr
private

Definition at line 167 of file fluid_structure_eigenproblem.cpp.

◆ matDFViscFluidPtr

boost::shared_ptr<MatrixDouble> FSI::matDFViscFluidPtr
private

Definition at line 168 of file fluid_structure_eigenproblem.cpp.

◆ matDFViscSolidPtr

boost::shared_ptr<MatrixDouble> FSI::matDFViscSolidPtr
private

Definition at line 169 of file fluid_structure_eigenproblem.cpp.

◆ matDSolidPtr

boost::shared_ptr<MatrixDouble> FSI::matDSolidPtr
private

Definition at line 166 of file fluid_structure_eigenproblem.cpp.

◆ matGradPtr

boost::shared_ptr<MatrixDouble> FSI::matGradPtr
private

Definition at line 163 of file fluid_structure_eigenproblem.cpp.

◆ matStrainPtr

boost::shared_ptr<MatrixDouble> FSI::matStrainPtr
private

Definition at line 164 of file fluid_structure_eigenproblem.cpp.

◆ matStressPtr

boost::shared_ptr<MatrixDouble> FSI::matStressPtr
private

Definition at line 165 of file fluid_structure_eigenproblem.cpp.

◆ mField

MoFEM::Interface& FSI::mField
private

Definition at line 153 of file fluid_structure_eigenproblem.cpp.

◆ pEP

SmartPetscObj<PEP> FSI::pEP
private

Definition at line 176 of file fluid_structure_eigenproblem.cpp.

◆ rFluid

Range FSI::rFluid
private

Definition at line 181 of file fluid_structure_eigenproblem.cpp.

◆ rInterface

Range FSI::rInterface
private

Definition at line 181 of file fluid_structure_eigenproblem.cpp.

◆ rSolid

Range FSI::rSolid
private

Definition at line 181 of file fluid_structure_eigenproblem.cpp.

◆ solverType

int FSI::solverType
private

Definition at line 179 of file fluid_structure_eigenproblem.cpp.


The documentation for this struct was generated from the following file: