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;
};
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 lo_uid =
DofEntity::getUniqueIdCalculate(0, get_id_for_min_type<MBVERTEX>());
auto hi_uid =
DofEntity::getUniqueIdCalculate(2, get_id_for_max_type<MBVERTEX>());
auto hi = dofs->upper_bound(lo_uid);
std::array<double, 3> coords;
for (auto lo = dofs->lower_bound(lo_uid); lo != hi; ++lo) {
auto ent = (*lo)->getEnt();
if ((*lo)->getDofCoeffIdx() == 0) {
INSERT_VALUES);
-coords[1], INSERT_VALUES);
-coords[2], INSERT_VALUES);
} else if ((*lo)->getDofCoeffIdx() == 1) {
INSERT_VALUES);
coords[0], INSERT_VALUES);
-coords[2], INSERT_VALUES);
} else if ((*lo)->getDofCoeffIdx() == 2) {
1, INSERT_VALUES);
coords[0], INSERT_VALUES);
coords[1], INSERT_VALUES);
}
}
}
}
SCATTER_FORWARD);
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(
};
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) {
}
CHKERR EPSSetDeflationSpace(
ePS, 6, &deflate_vectors[0]);
};
auto print_info = [&]() {
ST st;
PetscInt nev, maxit, its;
" Number of iterations of the method: %d", its);
MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
" Solution method: %s", type);
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 EPSSetWhichEigenpairs(
ePS, EPS_SMALLEST_MAGNITUDE);
};
}
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(),
OpPPMap::DataMapVec{},
OpPPMap::DataMapMat{{"U", u_ptr}},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{{"STRAIN", strain_ptr}, {"STRESS", stress_ptr}}
)
);
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[]) {
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
SlepcFinalize();
}
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
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.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'n', SPACE_DIM > n
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.
auto createDMVector(DM dm)
Get smart vector from DM.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
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)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double young_modulus
Young modulus.
#define EXECUTABLE_DIMENSION
static constexpr int approx_order
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Set up problem]
boost::shared_ptr< MatrixDouble > matDPtr
std::array< SmartPetscObj< Vec >, 6 > rigidBodyMotion
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode outputResults()
[Solve]
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
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.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Section manager is used to create indexes and sections.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
Set inverse jacobian to base functions.
PipelineManager interface.
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.