Calculate natural frequencies in 2d and 3d problems.
#undef EPS
#include <slepceps.h>
using DomainEle = FaceElementForcesAndSourcesCoreBase;
};
};
GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
private:
boost::shared_ptr<MatrixDouble> matGradPtr;
boost::shared_ptr<MatrixDouble> matStrainPtr;
boost::shared_ptr<MatrixDouble> matStressPtr;
boost::shared_ptr<MatrixDouble> matDPtr;
SmartPetscObj<EPS> ePS;
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;
};
matGradPtr = boost::make_shared<MatrixDouble>();
matStrainPtr = boost::make_shared<MatrixDouble>();
matStressPtr = boost::make_shared<MatrixDouble>();
matDPtr = boost::make_shared<MatrixDouble>();
matDPtr->resize(size_symm * size_symm, 1);
CHKERR set_matrial_stiffens();
}
}
auto simple = mField.getInterface<Simple>();
}
auto *
simple = mField.getInterface<Simple>();
}
auto *
simple = mField.getInterface<Simple>();
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 *
simple = mField.getInterface<Simple>();
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto calculate_stiffness_matrix = [&]() {
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetInvJacH1ForFace(
invJac));
}
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpK(
"U",
"U", matDPtr));
};
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
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 get_rho = [](
const double,
const double,
const double) {
return rho; };
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpMass(
"U",
"U", get_rho));
};
CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
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 simple = mField.getInterface<Simple>();
auto is_manager = mField.getInterface<ISManager>();
auto create_eps = [](MPI_Comm comm) {
return SmartPetscObj<EPS>(
eps);
};
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());
}
auto *pipeline_mng = mField.getInterface<PipelineManager>();
auto *
simple = mField.getInterface<Simple>();
pipeline_mng->getDomainLhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(
invJac));
}
post_proc_fe->getOpPtrVector().push_back(
matGradPtr));
post_proc_fe->getOpPtrVector().push_back(
new OpSymmetrizeTensor<SPACE_DIM>("U", matGradPtr, matStrainPtr));
post_proc_fe->getOpPtrVector().push_back(
new OpTensorTimesSymmetricTensor<SPACE_DIM, SPACE_DIM>(
"U", matStrainPtr, matStressPtr, matDPtr));
post_proc_fe->getOpPtrVector().push_back(new OpPostProcElastic<SPACE_DIM>(
"U", post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts, matStrainPtr,
matStressPtr));
post_proc_fe->addFieldValuesPostProc("U");
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, sqrt(std::abs(eigr)), 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[]) {
auto core_log = logging::core::get();
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
SlepcFinalize();
}
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Kronecker Delta class symmetric.
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
int main(int argc, char *argv[])
static char help[]
[Check]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
DataForcesAndSourcesCore::EntData EntData
auto smartCreateDMVector
Get smart vector from DM.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMass
VolumeElementForcesAndSourcesCore DomainEle
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
SmartPetscObj< Mat > smartMatDuplicate(Mat &mat, MatDuplicateOption op)
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
#define EXECUTABLE_DIMENSION
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
static constexpr int approx_order
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Applying essential BC]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode runProblem()
[Run problem]
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode outputResults()
[Solve]
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Data on single entity (This is passed as argument to DataOperator::doWork)
Deprecated interface functions.
static UId getUniqueIdCalculate(const DofIdx dof, UId ent_uid)
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
default operator for TET element