Calculate natural frequencies in 2d and 3d problems.
#undef EPS
#include <slepceps.h>
};
};
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
private:
boost::shared_ptr<MatrixDouble> matDPtr;
std::array<SmartPetscObj<Vec>, 6> rigidBodyMotion;
};
PETSC_NULL);
PETSC_NULL);
auto set_matrial_stiffens = [&]() {
auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
: 1;
};
matDPtr = boost::make_shared<MatrixDouble>();
CHKERR set_matrial_stiffens();
}
}
}
}
for (
int n = 1;
n != 6; ++
n)
auto problem_ptr = mField.get_problem(
simple->getProblemName());
auto dofs = problem_ptr->getNumeredRowDofsPtr();
const auto bit_number = mField.get_field_bit_number("U");
auto lo_uid =
auto hi_uid =
auto hi = dofs->upper_bound(lo_uid);
std::array<double, 3> coords;
for (auto lo = dofs->lower_bound(lo_uid); lo != hi; ++lo) {
if ((*lo)->getPart() == mField.get_comm_rank()) {
auto ent = (*lo)->getEnt();
CHKERR mField.get_moab().get_coords(&ent, 1, coords.data());
if ((*lo)->getDofCoeffIdx() == 0) {
CHKERR VecSetValue(rigidBodyMotion[0], (*lo)->getPetscGlobalDofIdx(), 1,
INSERT_VALUES);
CHKERR VecSetValue(rigidBodyMotion[3], (*lo)->getPetscGlobalDofIdx(),
-coords[1], INSERT_VALUES);
CHKERR VecSetValue(rigidBodyMotion[4], (*lo)->getPetscGlobalDofIdx(),
-coords[2], INSERT_VALUES);
} else if ((*lo)->getDofCoeffIdx() == 1) {
CHKERR VecSetValue(rigidBodyMotion[1], (*lo)->getPetscGlobalDofIdx(), 1,
INSERT_VALUES);
CHKERR VecSetValue(rigidBodyMotion[3], (*lo)->getPetscGlobalDofIdx(),
coords[0], INSERT_VALUES);
CHKERR VecSetValue(rigidBodyMotion[5], (*lo)->getPetscGlobalDofIdx(),
-coords[2], INSERT_VALUES);
} else if ((*lo)->getDofCoeffIdx() == 2) {
CHKERR VecSetValue(rigidBodyMotion[2], (*lo)->getPetscGlobalDofIdx(),
1, INSERT_VALUES);
CHKERR VecSetValue(rigidBodyMotion[4], (*lo)->getPetscGlobalDofIdx(),
coords[0], INSERT_VALUES);
CHKERR VecSetValue(rigidBodyMotion[5], (*lo)->getPetscGlobalDofIdx(),
coords[1], INSERT_VALUES);
}
}
}
}
for (
int n = 0;
n != rigidBodyMotion.size(); ++
n) {
CHKERR VecAssemblyBegin(rigidBodyMotion[
n]);
CHKERR VecAssemblyEnd(rigidBodyMotion[
n]);
CHKERR VecGhostUpdateBegin(rigidBodyMotion[
n], INSERT_VALUES,
SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(rigidBodyMotion[
n], INSERT_VALUES,
SCATTER_FORWARD);
}
}
auto calculate_stiffness_matrix = [&]() {
pipeline_mng->getDomainLhsFE().reset();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpK(
"U",
"U", matDPtr));
};
pipeline_mng->getDomainLhsFE()->B =
K;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(
K, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
K, MAT_FINAL_ASSEMBLY);
};
auto calculate_mass_matrix = [&]() {
pipeline_mng->getDomainLhsFE().reset();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass(
"U",
"U", get_rho));
};
pipeline_mng->getDomainLhsFE()->B =
M;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
};
CHKERR calculate_stiffness_matrix();
CHKERR calculate_mass_matrix();
}
auto create_eps = [](MPI_Comm comm) {
};
auto deflate_vectors = [&]() {
std::array<Vec, 6> deflate_vectors;
for (
int n = 0;
n != 6; ++
n) {
deflate_vectors[
n] = rigidBodyMotion[
n];
}
CHKERR EPSSetDeflationSpace(ePS, 6, &deflate_vectors[0]);
};
auto print_info = [&]() {
ST st;
PetscInt nev, maxit, its;
CHKERR EPSGetIterationNumber(ePS, &its);
" Number of iterations of the method: %d", its);
CHKERR EPSGetDimensions(ePS, &nev, NULL, NULL);
MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
" Number of requested eigenvalues: %d",
nev);
" Stopping condition: tol=%.4g, maxit=%d", (
double)
tol, maxit);
PetscScalar eigr, eigi;
for (int nn = 0; nn < nev; nn++) {
CHKERR EPSGetEigenpair(ePS, nn, &eigr, &eigi, PETSC_NULL, PETSC_NULL);
" ncov = %d eigr = %.4g eigi = %.4g (inv eigr = %.4g)", nn,
eigr, eigi, 1. / eigr);
}
};
auto setup_eps = [&]() {
CHKERR EPSSetProblemType(ePS, EPS_GHEP);
CHKERR EPSSetWhichEigenpairs(ePS, EPS_SMALLEST_MAGNITUDE);
CHKERR EPSSetFromOptions(ePS);
};
ePS = create_eps(mField.get_comm());
}
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<MatrixDouble>();
auto grad_ptr = boost::make_shared<MatrixDouble>();
auto strain_ptr = boost::make_shared<MatrixDouble>();
auto stress_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
"U", strain_ptr, stress_ptr, matDPtr));
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
)
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
PetscInt nev;
CHKERR EPSGetDimensions(ePS, &nev, NULL, NULL);
PetscScalar eigr, eigi, nrm2r;
for (int nn = 0; nn < nev; nn++) {
CHKERR EPSGetEigenpair(ePS, nn, &eigr, &eigi,
D, PETSC_NULL);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
" ncov = %d omega2 = %.8g omega = %.8g frequency = %.8g", nn,
eigr, std::sqrt(std::abs(eigr)),
std::sqrt(std::abs(eigr)) / (2 * M_PI));
CHKERR pipeline_mng->loopFiniteElements();
post_proc_fe->writeFile("out_eig_" + boost::lexical_cast<std::string>(nn) +
".h5m");
}
}
PetscBool test_flg = PETSC_FALSE;
if (test_flg) {
PetscScalar eigr, eigi;
CHKERR EPSGetEigenpair(
ePS, 0, &eigr, &eigi, PETSC_NULL, PETSC_NULL);
constexpr double regression_value = 12579658;
if (fabs(eigr - regression_value) > 1)
"Regression test faileed; wrong eigen value.");
}
}
static char help[] =
"...\n\n";
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
SlepcInitialize(&argc, &argv, param_file,
help);
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
}
SlepcFinalize();
}