692 {
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(
707 post_proc_fe->getOpPtrVector().push_back(
709 post_proc_fe->getOpPtrVector().push_back(
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
719
720 post_proc_fe->getOpPtrVector().push_back(
722 post_proc_fe->getOpPtrVector().push_back(
724 post_proc_fe->getOpPtrVector().push_back(
726
727 post_proc_fe->getOpPtrVector().push_back(
729 post_proc_fe->getOpPtrVector().push_back(
731 post_proc_fe->getOpPtrVector().push_back(
733
735
736 post_proc_fe->getOpPtrVector().push_back(
737
738 new OpPPMap(post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts,
739
741
743
745
747 {"STRAINf", strain_f_ptr}}
748
749 )
750
751 );
752
755
757 op_tag_handle->doWorkRhsHook = [&](
DataOperator *op_ptr,
int side,
761
762 auto get_tag = [&]() {
763 int def_val = -1;
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");
770 };
771
772 auto fe_ent =
static_cast<DomainEleOp *
>(op_ptr)->getFEEntityHandle();
773
777 }
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) {
800 CHKERR pipeline_mng->loopFiniteElements();
801 post_proc_fe->writeFile(file_name);
803 };
804
805 auto post_proc_eps = [&](
auto eps) {
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
824 "out_eig_" + boost::lexical_cast<std::string>(nn) +
".h5m",
v);
825 }
826
828 };
829
830 auto post_proc_pep = [&](
auto eps) {
832
835
836 PetscInt nev;
838 PetscScalar ki, kr, error, nrm2r;
839
840 int idx = 0;
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
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
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
875
877}
#define MOFEM_LOG_C(channel, severity, format,...)
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
ElementsAndOps< SPACE_DIM >::DomainEleOp DomainEleOp
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
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)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
BcMapByBlockName & getBcMapByBlockName()
Get the bc map.
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.
