34 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
84 auto set_matrial_stiffens = [&]() {
91 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*
matDPtr);
104 matDPtr = boost::make_shared<MatrixDouble>();
109 CHKERR set_matrial_stiffens();
163 for (
int n = 1;
n != 6; ++
n)
177 auto hi = dofs->upper_bound(lo_uid);
178 std::array<double, 3> coords;
180 for (
auto lo = dofs->lower_bound(lo_uid); lo != hi; ++lo) {
184 auto ent = (*lo)->getEnt();
187 if ((*lo)->getDofCoeffIdx() == 0) {
191 -coords[1], INSERT_VALUES);
194 -coords[2], INSERT_VALUES);
196 }
else if ((*lo)->getDofCoeffIdx() == 1) {
200 coords[0], INSERT_VALUES);
203 -coords[2], INSERT_VALUES);
205 }
else if ((*lo)->getDofCoeffIdx() == 2) {
210 coords[0], INSERT_VALUES);
212 coords[1], INSERT_VALUES);
237 auto dm =
simple->getDM();
241 auto calculate_stiffness_matrix = [&]() {
243 pipeline_mng->getDomainLhsFE().reset();
245 auto det_ptr = boost::make_shared<VectorDouble>();
246 auto jac_ptr = boost::make_shared<MatrixDouble>();
247 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
248 pipeline_mng->getOpDomainLhsPipeline().push_back(
250 pipeline_mng->getOpDomainLhsPipeline().push_back(
252 pipeline_mng->getOpDomainLhsPipeline().push_back(
254 pipeline_mng->getOpDomainLhsPipeline().push_back(
257 pipeline_mng->getOpDomainLhsPipeline().push_back(
264 pipeline_mng->getDomainLhsFE()->B =
K;
266 CHKERR pipeline_mng->loopFiniteElements();
267 CHKERR MatAssemblyBegin(
K, MAT_FINAL_ASSEMBLY);
268 CHKERR MatAssemblyEnd(
K, MAT_FINAL_ASSEMBLY);
272 auto calculate_mass_matrix = [&]() {
274 pipeline_mng->getDomainLhsFE().reset();
276 auto det_ptr = boost::make_shared<VectorDouble>();
277 auto jac_ptr = boost::make_shared<MatrixDouble>();
278 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
279 pipeline_mng->getOpDomainLhsPipeline().push_back(
281 pipeline_mng->getOpDomainLhsPipeline().push_back(
283 pipeline_mng->getOpDomainLhsPipeline().push_back(
287 pipeline_mng->getOpDomainLhsPipeline().push_back(
288 new OpMass(
"U",
"U", get_rho));
294 pipeline_mng->getDomainLhsFE()->B =
M;
295 CHKERR pipeline_mng->loopFiniteElements();
296 CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
297 CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
301 CHKERR calculate_stiffness_matrix();
302 CHKERR calculate_mass_matrix();
314 auto create_eps = [](MPI_Comm comm) {
320 auto deflate_vectors = [&]() {
323 std::array<Vec, 6> deflate_vectors;
324 for (
int n = 0;
n != 6; ++
n) {
327 CHKERR EPSSetDeflationSpace(
ePS, 6, &deflate_vectors[0]);
331 auto print_info = [&]() {
336 PetscInt nev, maxit, its;
340 " Number of iterations of the method: %d", its);
343 MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
" Solution method: %s", type);
344 CHKERR EPSGetDimensions(
ePS, &nev, NULL, NULL);
345 MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
" Number of requested eigenvalues: %d",
349 " Stopping condition: tol=%.4g, maxit=%d", (
double)
tol, maxit);
351 PetscScalar eigr, eigi;
352 for (
int nn = 0; nn < nev; nn++) {
353 CHKERR EPSGetEigenpair(
ePS, nn, &eigr, &eigi, PETSC_NULL, PETSC_NULL);
355 " ncov = %d eigr = %.4g eigi = %.4g (inv eigr = %.4g)", nn,
356 eigr, eigi, 1. / eigr);
362 auto setup_eps = [&]() {
365 CHKERR EPSSetWhichEigenpairs(
ePS, EPS_SMALLEST_MAGNITUDE);
396 pipeline_mng->getDomainLhsFE().reset();
397 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
399 auto det_ptr = boost::make_shared<VectorDouble>();
400 auto jac_ptr = boost::make_shared<MatrixDouble>();
401 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
402 post_proc_fe->getOpPtrVector().push_back(
404 post_proc_fe->getOpPtrVector().push_back(
406 post_proc_fe->getOpPtrVector().push_back(
409 auto u_ptr = boost::make_shared<MatrixDouble>();
410 auto grad_ptr = boost::make_shared<MatrixDouble>();
411 auto strain_ptr = boost::make_shared<MatrixDouble>();
412 auto stress_ptr = boost::make_shared<MatrixDouble>();
414 post_proc_fe->getOpPtrVector().push_back(
416 post_proc_fe->getOpPtrVector().push_back(
418 post_proc_fe->getOpPtrVector().push_back(
420 post_proc_fe->getOpPtrVector().push_back(
422 "U", strain_ptr, stress_ptr,
matDPtr));
426 post_proc_fe->getOpPtrVector().push_back(
429 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
443 pipeline_mng->getDomainRhsFE() = post_proc_fe;
445 auto dm =
simple->getDM();
449 CHKERR EPSGetDimensions(ePS, &nev, NULL, NULL);
450 PetscScalar eigr, eigi, nrm2r;
451 for (
int nn = 0; nn < nev; nn++) {
452 CHKERR EPSGetEigenpair(ePS, nn, &eigr, &eigi,
D, PETSC_NULL);
453 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
454 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
455 CHKERR VecNorm(
D, NORM_2, &nrm2r);
457 " ncov = %d omega2 = %.8g omega = %.8g frequency = %.8g", nn,
458 eigr, std::sqrt(std::abs(eigr)),
459 std::sqrt(std::abs(eigr)) / (2 * M_PI));
461 CHKERR pipeline_mng->loopFiniteElements();
462 post_proc_fe->writeFile(
"out_eig_" + boost::lexical_cast<std::string>(nn) +
473 PetscBool test_flg = PETSC_FALSE;
476 PetscScalar eigr, eigi;
477 CHKERR EPSGetEigenpair(
ePS, 0, &eigr, &eigi, PETSC_NULL, PETSC_NULL);
478 constexpr double regression_value = 12579658;
479 if (fabs(eigr - regression_value) > 1)
481 "Regression test faileed; wrong eigen value.");
489int main(
int argc,
char *argv[]) {
497 auto core_log = logging::core::get();
506 DMType dm_name =
"DMMOFEM";
511 moab::Core mb_instance;
512 moab::Interface &moab = mb_instance;
#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.
static char help[]
[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.
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< '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()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
boost::shared_ptr< MatrixDouble > matDPtr
std::array< SmartPetscObj< Vec >, 6 > rigidBodyMotion
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
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.
static UId getUniqueIdCalculate(const DofIdx dof, UId ent_uid)
Data on single entity (This is passed as argument to DataOperator::doWork)
Section manager is used to create indexes and sections.
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.
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
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.