25 using namespace MoFEM;
57 GAUSS>::OpGradSymTensorGrad<1, SPACE_DIM, SPACE_DIM, 0>;
59 GAUSS>::OpBaseTimesVector<1, SPACE_DIM, 0>;
89 boost::shared_ptr<MatrixDouble> matGradPtr;
90 boost::shared_ptr<MatrixDouble> matStrainPtr;
91 boost::shared_ptr<MatrixDouble> matStressPtr;
92 boost::shared_ptr<MatrixDouble> matDPtr;
101 auto set_material_stiffness = [&]() {
112 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*matDPtr);
121 auto set_body_force = [&]() {
124 auto t_force = getFTensor1FromMat<SPACE_DIM, 0>(*bodyForceMatPtr);
132 matGradPtr = boost::make_shared<MatrixDouble>();
133 matStrainPtr = boost::make_shared<MatrixDouble>();
134 matStressPtr = boost::make_shared<MatrixDouble>();
135 matDPtr = boost::make_shared<MatrixDouble>();
136 bodyForceMatPtr = boost::make_shared<MatrixDouble>();
139 matDPtr->resize(size_symm * size_symm, 1);
143 CHKERR set_material_stiffness();
155 CHKERR createCommonData();
156 CHKERR boundaryCondition();
197 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_X",
199 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Y",
201 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_Z",
203 CHKERR bc_mng->removeBlockDOFsOnEntities(
simple->getProblemName(),
"FIX_ALL",
207 pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBodyForce(
208 "U", bodyForceMatPtr, [](
double,
double,
double) {
return 1.; }));
221 auto det_ptr = boost::make_shared<VectorDouble>();
222 auto jac_ptr = boost::make_shared<MatrixDouble>();
223 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
248 auto solver = pipeline_mng->createKSP();
249 CHKERR KSPSetFromOptions(solver);
252 auto dm =
simple->getDM();
257 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
258 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
268 auto det_ptr = boost::make_shared<VectorDouble>();
269 auto jac_ptr = boost::make_shared<MatrixDouble>();
270 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
272 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
273 post_proc_fe->generateReferenceElementMesh();
275 post_proc_fe->getOpPtrVector().push_back(
277 post_proc_fe->getOpPtrVector().push_back(
279 post_proc_fe->getOpPtrVector().push_back(
282 post_proc_fe->getOpPtrVector().push_back(
285 post_proc_fe->getOpPtrVector().push_back(
287 post_proc_fe->getOpPtrVector().push_back(
289 "U", matStrainPtr, matStressPtr, matDPtr));
291 "U", post_proc_fe->postProcMesh, post_proc_fe->mapGaussPts, matStrainPtr,
293 post_proc_fe->addFieldValuesPostProc(
"U");
296 CHKERR post_proc_fe->writeFile(
"out_elastic.h5m");
307 pipeline_mng->getDomainRhsFE().reset();
308 pipeline_mng->getDomainLhsFE().reset();
311 auto det_ptr = boost::make_shared<VectorDouble>();
312 auto jac_ptr = boost::make_shared<MatrixDouble>();
313 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
314 pipeline_mng->getOpDomainRhsPipeline().push_back(
316 pipeline_mng->getOpDomainRhsPipeline().push_back(
318 pipeline_mng->getOpDomainRhsPipeline().push_back(
321 pipeline_mng->getOpDomainRhsPipeline().push_back(
324 pipeline_mng->getOpDomainRhsPipeline().push_back(
326 pipeline_mng->getOpDomainRhsPipeline().push_back(
328 "U", matStrainPtr, matStressPtr, matDPtr));
329 pipeline_mng->getOpDomainRhsPipeline().push_back(
331 pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpBodyForce(
332 "U", bodyForceMatPtr, [](
double,
double,
double) {
return -1.; }));
337 auto dm =
simple->getDM();
339 pipeline_mng->getDomainRhsFE()->ksp_f = res;
341 CHKERR VecZeroEntries(res);
342 CHKERR pipeline_mng->loopFiniteElements();
343 CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
344 CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
345 CHKERR VecAssemblyBegin(res);
346 CHKERR VecAssemblyEnd(res);
349 CHKERR VecNorm(res, NORM_2, &nrm2);
350 MOFEM_LOG_C(
"WORLD", Sev::verbose,
"residual = %3.4e\n", nrm2);
351 constexpr
double eps = 1e-8;
359 static char help[] =
"...\n\n";
361 int main(
int argc,
char *argv[]) {
370 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)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Kronecker Delta class symmetric.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForce
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, SPACE_DIM > OpBodyForce
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#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[])
EntitiesFieldData::EntData EntData
[Define dimension]
static char help[]
[Check]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 0 > OpBodyForce
constexpr int SPACE_DIM
[Define dimension]
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
constexpr double young_modulus
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 DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset 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)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
static constexpr int approx_order
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
PipelineManager::FaceEle DomainEle
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
boost::shared_ptr< MatrixDouble > bodyForceMatPtr
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
Simple interface for fast problem set-up.
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.
default operator for EDGE element
Data on single entity (This is passed as argument to DataOperator::doWork)
default operator for TRI element
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element with switches.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator