34 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
84 auto set_matrial_stiffens = [&]() {
104 matDPtr = boost::make_shared<MatrixDouble>();
109 CHKERR set_matrial_stiffens();
163 for (
int n = 1;
n != 6; ++
n)
174 auto hi = dofs->upper_bound(lo_uid);
175 std::array<double, 3> coords;
177 for (
auto lo = dofs->lower_bound(lo_uid); lo != hi; ++lo) {
181 auto ent = (*lo)->getEnt();
184 if ((*lo)->getDofCoeffIdx() == 0) {
188 -coords[1], INSERT_VALUES);
191 -coords[2], INSERT_VALUES);
193 }
else if ((*lo)->getDofCoeffIdx() == 1) {
197 coords[0], INSERT_VALUES);
200 -coords[2], INSERT_VALUES);
202 }
else if ((*lo)->getDofCoeffIdx() == 2) {
207 coords[0], INSERT_VALUES);
209 coords[1], INSERT_VALUES);
234 auto dm =
simple->getDM();
238 auto calculate_stiffness_matrix = [&]() {
240 pipeline_mng->getDomainLhsFE().reset();
242 auto det_ptr = boost::make_shared<VectorDouble>();
243 auto jac_ptr = boost::make_shared<MatrixDouble>();
244 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
245 pipeline_mng->getOpDomainLhsPipeline().push_back(
247 pipeline_mng->getOpDomainLhsPipeline().push_back(
249 pipeline_mng->getOpDomainLhsPipeline().push_back(
251 pipeline_mng->getOpDomainLhsPipeline().push_back(
254 pipeline_mng->getOpDomainLhsPipeline().push_back(
261 pipeline_mng->getDomainLhsFE()->B =
K;
263 CHKERR pipeline_mng->loopFiniteElements();
264 CHKERR MatAssemblyBegin(
K, MAT_FINAL_ASSEMBLY);
265 CHKERR MatAssemblyEnd(
K, MAT_FINAL_ASSEMBLY);
269 auto calculate_mass_matrix = [&]() {
271 pipeline_mng->getDomainLhsFE().reset();
273 auto det_ptr = boost::make_shared<VectorDouble>();
274 auto jac_ptr = boost::make_shared<MatrixDouble>();
275 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
276 pipeline_mng->getOpDomainLhsPipeline().push_back(
278 pipeline_mng->getOpDomainLhsPipeline().push_back(
280 pipeline_mng->getOpDomainLhsPipeline().push_back(
284 pipeline_mng->getOpDomainLhsPipeline().push_back(
285 new OpMass(
"U",
"U", get_rho));
291 pipeline_mng->getDomainLhsFE()->B =
M;
292 CHKERR pipeline_mng->loopFiniteElements();
293 CHKERR MatAssemblyBegin(
M, MAT_FINAL_ASSEMBLY);
294 CHKERR MatAssemblyEnd(
M, MAT_FINAL_ASSEMBLY);
298 CHKERR calculate_stiffness_matrix();
299 CHKERR calculate_mass_matrix();
309 auto create_eps = [](MPI_Comm comm) {
315 auto deflate_vectors = [&]() {
318 std::array<Vec, 6> deflate_vectors;
319 for (
int n = 0;
n != 6; ++
n) {
322 CHKERR EPSSetDeflationSpace(
ePS, 6, &deflate_vectors[0]);
326 auto print_info = [&]() {
331 PetscInt nev, maxit, its;
335 " Number of iterations of the method: %d", its);
338 MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
" Solution method: %s", type);
339 CHKERR EPSGetDimensions(
ePS, &nev, NULL, NULL);
340 MOFEM_LOG_C(
"EXAMPLE", Sev::inform,
" Number of requested eigenvalues: %d",
344 " Stopping condition: tol=%.4g, maxit=%d", (
double)
tol, maxit);
346 PetscScalar eigr, eigi;
347 for (
int nn = 0; nn < nev; nn++) {
348 CHKERR EPSGetEigenpair(
ePS, nn, &eigr, &eigi, PETSC_NULLPTR, PETSC_NULLPTR);
350 " ncov = %d eigr = %.4g eigi = %.4g (inv eigr = %.4g)", nn,
351 eigr, eigi, 1. / eigr);
357 auto setup_eps = [&]() {
360 CHKERR EPSSetWhichEigenpairs(
ePS, EPS_SMALLEST_MAGNITUDE);
391 pipeline_mng->getDomainLhsFE().reset();
392 auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
394 auto det_ptr = boost::make_shared<VectorDouble>();
395 auto jac_ptr = boost::make_shared<MatrixDouble>();
396 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
397 post_proc_fe->getOpPtrVector().push_back(
399 post_proc_fe->getOpPtrVector().push_back(
401 post_proc_fe->getOpPtrVector().push_back(
404 auto u_ptr = boost::make_shared<MatrixDouble>();
405 auto grad_ptr = boost::make_shared<MatrixDouble>();
406 auto strain_ptr = boost::make_shared<MatrixDouble>();
407 auto stress_ptr = boost::make_shared<MatrixDouble>();
409 post_proc_fe->getOpPtrVector().push_back(
411 post_proc_fe->getOpPtrVector().push_back(
413 post_proc_fe->getOpPtrVector().push_back(
415 post_proc_fe->getOpPtrVector().push_back(
417 strain_ptr, stress_ptr,
matDPtr));
421 post_proc_fe->getOpPtrVector().push_back(
424 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
438 pipeline_mng->getDomainRhsFE() = post_proc_fe;
440 auto dm =
simple->getDM();
444 CHKERR EPSGetDimensions(ePS, &nev, NULL, NULL);
445 PetscScalar eigr, eigi, nrm2r;
446 for (
int nn = 0; nn < nev; nn++) {
447 CHKERR EPSGetEigenpair(ePS, nn, &eigr, &eigi,
D, PETSC_NULLPTR);
448 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
449 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
450 CHKERR VecNorm(
D, NORM_2, &nrm2r);
452 " ncov = %d omega2 = %.8g omega = %.8g frequency = %.8g", nn,
453 eigr, std::sqrt(std::abs(eigr)),
454 std::sqrt(std::abs(eigr)) / (2 * M_PI));
456 CHKERR pipeline_mng->loopFiniteElements();
457 post_proc_fe->writeFile(
"out_eig_" + boost::lexical_cast<std::string>(nn) +
468 PetscBool test_flg = PETSC_FALSE;
471 PetscScalar eigr, eigi;
472 CHKERR EPSGetEigenpair(
ePS, 0, &eigr, &eigi, PETSC_NULLPTR, PETSC_NULLPTR);
473 constexpr double regression_value = 12579658;
474 if (fabs(eigr - regression_value) > 1)
476 "Regression test faileed; wrong eigen value.");
484int main(
int argc,
char *argv[]) {
487 const char param_file[] =
"param_file.petsc";
488 SlepcInitialize(&argc, &argv, param_file,
help);
492 auto core_log = logging::core::get();
501 DMType dm_name =
"DMMOFEM";
506 moab::Core mb_instance;
507 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)
constexpr int SPACE_DIM
[Define dimension]
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
constexpr double young_modulus
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
[Define dimension]
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 >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
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< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
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
EntityHandle get_id_for_min_type()
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.
static FTensor::Ddg< FTensor::PackPtr< T *, 1 >, Tensor_Dim01, Tensor_Dim23 > getFTensor4DdgFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 4 on first two and last indices from form data matrix.
SmartPetscObj< Mat > matDuplicate(Mat mat, MatDuplicateOption op)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
#define EXECUTABLE_DIMENSION
static constexpr int approx_order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
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 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)
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 field rank 0, i.e. vector field.
Get values at integration pts for tensor field rank 1, i.e. vector field.
Post post-proc data at points from hash maps.
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
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.
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Volume finite element base.