25 using namespace MoFEM;
49 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
55 GAUSS>::OpSource<1, SPACE_DIM>;
59 constexpr
double rho = 1;
91 boost::shared_ptr<MatrixDouble> matGradPtr;
92 boost::shared_ptr<MatrixDouble> matStrainPtr;
93 boost::shared_ptr<MatrixDouble> matStressPtr;
96 boost::shared_ptr<MatrixDouble> matDPtr;
103 auto set_matrial_stiffens = [&]() {
114 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
117 t_kd(
i,
j) * t_kd(
k,
l);
121 matGradPtr = boost::make_shared<MatrixDouble>();
122 matStrainPtr = boost::make_shared<MatrixDouble>();
123 matStressPtr = boost::make_shared<MatrixDouble>();
124 matAccelerationPtr = boost::make_shared<MatrixDouble>();
125 matInertiaPtr = boost::make_shared<MatrixDouble>();
126 matDPtr = boost::make_shared<MatrixDouble>();
129 matDPtr->resize(size_symm * size_symm, 1);
131 CHKERR set_matrial_stiffens();
142 CHKERR createCommonData();
143 CHKERR boundaryCondition();
181 auto fix_disp = [&](
const std::string blockset_name) {
184 if (it->getName().compare(0, blockset_name.length(), blockset_name) ==
186 CHKERR mField.get_moab().get_entities_by_handle(it->meshset, fix_ents,
193 auto remove_ents = [&](
const Range &&ents,
const int lo,
const int hi) {
198 CHKERR mField.get_moab().get_connectivity(ents, verts,
true);
202 CHKERR mField.get_moab().get_adjacencies(ents, 1,
false, adj,
203 moab::Interface::UNION);
207 CHKERR prb_mng->removeDofsOnEntities(
simple->getProblemName(),
"U", verts,
212 CHKERR remove_ents(fix_disp(
"FIX_X"), 0, 0);
213 CHKERR remove_ents(fix_disp(
"FIX_Y"), 1, 1);
214 CHKERR remove_ents(fix_disp(
"FIX_Z"), 2, 2);
215 CHKERR remove_ents(fix_disp(
"FIX_ALL"), 0, 3);
228 pipeline_mng->getOpDomainLhsPipeline().push_back(
230 pipeline_mng->getOpDomainLhsPipeline().push_back(
232 pipeline_mng->getOpDomainRhsPipeline().push_back(
234 pipeline_mng->getOpDomainRhsPipeline().push_back(
239 auto get_rho = [
this](
const double,
const double,
const double) {
241 auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
242 return rho * fe_domain_lhs->ts_aa;
245 auto get_body_force = [
this](
const double,
const double,
const double) {
247 auto fe_domain_rhs = pipeline_mng->getDomainRhsFE();
250 const auto time = fe_domain_rhs->ts_t;
251 t_source(
i) *= sin(time *
omega * M_PI);
255 pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpK(
"U",
"U", matDPtr));
256 pipeline_mng->getOpDomainLhsPipeline().push_back(
257 new OpMass(
"U",
"U", get_rho));
259 pipeline_mng->getOpDomainRhsPipeline().push_back(
261 pipeline_mng->getOpDomainRhsPipeline().push_back(
264 pipeline_mng->getOpDomainRhsPipeline().push_back(
266 pipeline_mng->getOpDomainRhsPipeline().push_back(
268 "U", matStrainPtr, matStressPtr, matDPtr));
269 pipeline_mng->getOpDomainRhsPipeline().push_back(
271 pipeline_mng->getOpDomainRhsPipeline().push_back(
273 matAccelerationPtr));
274 pipeline_mng->getOpDomainRhsPipeline().push_back(
276 pipeline_mng->getOpDomainRhsPipeline().push_back(
282 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
283 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
297 : dM(dm), postProc(post_proc){};
303 CHKERR postProc->writeFile(
304 "out_step_" + boost::lexical_cast<std::string>(ts_step) +
".h5m");
320 auto dm =
simple->getDM();
321 auto ts = pipeline_mng->createTS2();
324 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
325 post_proc_fe->generateReferenceElementMesh();
327 post_proc_fe->getOpPtrVector().push_back(
331 post_proc_fe->getOpPtrVector().push_back(
334 post_proc_fe->getOpPtrVector().push_back(
336 post_proc_fe->getOpPtrVector().push_back(
338 "U", matStrainPtr, matStressPtr, matDPtr));
340 "U", post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts, matStrainPtr,
342 post_proc_fe->addFieldValuesPostProc(
"U");
345 boost::shared_ptr<FEMethod> null_fe;
346 auto monitor_ptr = boost::make_shared<Monitor>(dm, post_proc_fe);
348 null_fe, monitor_ptr);
351 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
352 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
359 CHKERR TS2SetSolution(ts,
T, TT);
360 CHKERR TSSetFromOptions(ts);
363 CHKERR TSGetTime(ts, &ftime);
365 PetscInt steps, snesfails, rejects, nonlinits, linits;
366 CHKERR TSGetTimeStepNumber(ts, &steps);
367 CHKERR TSGetSNESFailures(ts, &snesfails);
368 CHKERR TSGetStepRejections(ts, &rejects);
369 CHKERR TSGetSNESIterations(ts, &nonlinits);
370 CHKERR TSGetKSPIterations(ts, &linits);
372 "steps %D (%D rejected, %D SNES fails), ftime %g, nonlinits "
374 steps, rejects, snesfails, ftime, nonlinits, linits);
383 PetscBool test_flg = PETSC_FALSE;
391 CHKERR VecNorm(
T, NORM_2, &nrm2);
392 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Regression norm " << nrm2;
393 constexpr
double regression_value = 1.09572;
394 if (fabs(nrm2 - regression_value) > 1e-2)
396 "Regression test faileed; wrong norm value.");
409 static char help[] =
"...\n\n";
411 int main(
int argc,
char *argv[]) {
418 auto core_log = logging::core::get();
427 DMType dm_name =
"DMMOFEM";
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
static double * ts_aa_ptr
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpK
constexpr double poisson_ratio
constexpr double shear_modulus_G
constexpr double bulk_modulus_K
DataForcesAndSourcesCore::EntData EntData
static double * ts_time_ptr
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
constexpr double young_modulus
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
auto smartCreateDMVector
Get smart vector from DM.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
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.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, 3, 1 > OpInertiaForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 3 > OpBodyForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 3 > OpMass
VolumeElementForcesAndSourcesCore DomainEle
const int save_every_nth_step
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
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)
OpCalculateInvJacForFaceImpl< 2 > OpCalculateInvJacForFace
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
ElementsAndOps< SPACE_DIM >::PostProcEle PostProcEle
static constexpr int approx_order
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
boost::shared_ptr< MatrixDouble > matAccelerationPtr
boost::shared_ptr< MatrixDouble > matInertiaPtr
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
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.
structure for User Loop Methods on finite elements
default operator for TRI element
Face finite element switched.
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.
Approximate field valuse for given petsc vector.
PipelineManager interface.
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
default operator for TET element
Volume finite element base.
Volume finite element with switches.
[Push operators to pipeline]
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcEle > postProc