#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();
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(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
}
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 get_rho = [](
const double,
const double,
const double) {
return rho; };
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 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();
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(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpSetInvJacH1ForFace(inv_jac_ptr));
}
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, 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(
try {
DMType dm_name = "DMMOFEM";
}
SlepcFinalize();
}
#define MOFEM_LOG_C(channel, severity, format,...)
void simple(double P1[], double P2[], double P3[], double c[], const int N)
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpK
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
constexpr double young_modulus
#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.
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< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
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)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
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
static constexpr int approx_order
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
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.
Deprecated interface functions.
static UId getUniqueIdCalculate(const DofIdx dof, UId ent_uid)
Data on single entity (This is passed as argument to DataOperator::doWork)
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.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator