v0.15.0
Loading...
Searching...
No Matches
Example Struct Reference

[Example] More...

Collaboration diagram for Example:
[legend]

Classes

struct  BoundaryOp
 
struct  CommonData
 [Example] More...
 
struct  DynamicFirstOrderConsConstantTimeScale
 
struct  DynamicFirstOrderConsSinusTimeScale
 
struct  OpCalcSurfaceAverageTemperature
 
struct  OpError
 
struct  OpError< 1 >
 
struct  OpFirst
 
struct  OpFluxRhs
 
struct  OpRadiationLhs
 
struct  OpRadiationRhs
 
struct  OpRhs
 
struct  OpSecond
 [Operator] More...
 
struct  OpZero
 [Common data] More...
 
struct  ScaledTimeScale
 

Public Types

enum  { VOL , COUNT }
 

Public Member Functions

 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 [Run problem]
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 
 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 

Static Public Attributes

static std::array< double, 2 > meshVolumeAndCount = {0, 0}
 

Private Types

enum  BoundingBox {
  CENTER_X = 0 , CENTER_Y , MAX_X , MAX_Y ,
  MIN_X , MIN_Y , LAST_BB
}
 

Private Member Functions

MoFEMErrorCode setupProblem ()
 [Run problem]
 
MoFEMErrorCode createCommonData ()
 [Set up problem]
 
MoFEMErrorCode bC ()
 [Create common data]
 
MoFEMErrorCode OPs ()
 [Boundary condition]
 
MoFEMErrorCode tsSolve ()
 
MoFEMErrorCode testOperators ()
 [Solve]
 
MoFEMErrorCode readMesh ()
 [Run problem]
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode boundaryCondition ()
 [Set up problem]
 
MoFEMErrorCode assembleSystem ()
 [Push operators to pipeline]
 
MoFEMErrorCode solveSystem ()
 [Solve]
 
MoFEMErrorCode outputResults ()
 [Solve]
 
MoFEMErrorCode checkResults ()
 [Postprocess results]
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode setUp ()
 [Run all]
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode setFieldValues ()
 [Create common data]
 
MoFEMErrorCode pushOperators ()
 [Set density distribution]
 
MoFEMErrorCode integrateElements ()
 [Push operators to pipeline]
 
MoFEMErrorCode postProcess ()
 [Integrate]
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 [Set up problem]
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode assembleSystemIntensity ()
 [Calculate flux on boundary]
 
MoFEMErrorCode assembleSystemFlux ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode calculateFlux (double &calc_flux)
 [Set up problem]
 
MoFEMErrorCode outputResults (const int i)
 [Solve]
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode bC ()
 
MoFEMErrorCode OPs ()
 
MoFEMErrorCode kspSolve ()
 [Push operators to pipeline]
 
MoFEMErrorCode postProcess ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setupAdJoint ()
 [Set up problem]
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode topologyModes ()
 [Boundary condition]
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveElastic ()
 [Push operators to pipeline]
 
MoFEMErrorCode postprocessElastic (int iter, SmartPetscObj< Vec > adjoint_vector=nullptr)
 [Solve]
 
MoFEMErrorCode calculateGradient (PetscReal *objective_function_value, Vec objective_function_gradient, Vec adjoint_vector)
 

Static Private Member Functions

static std::pair< int, int > getCoordsInImage (double x, double y)
 
static double rhsSource (const double x, const double y, const double)
 
static double lhsFlux (const double x, const double y, const double)
 
static int integrationRule (int, int, int p_data)
 

Private Attributes

MoFEM::InterfacemField
 
boost::shared_ptr< DomainElereactionFe
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
 
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
 
SimplesimpleInterface
 
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
 
boost::shared_ptr< CommonDatacommonDataPtr
 
FieldApproximationBase base
 
FieldSpace space
 
boost::shared_ptr< VectorDoubleapproxVals
 
boost::shared_ptr< MatrixDoubleapproxGradVals
 
Range pinchNodes
 
boost::shared_ptr< MatrixDoublevectorFieldPtr = nullptr
 
boost::shared_ptr< MatrixDoublematDPtr
 
SmartPetscObj< Mat > M
 
SmartPetscObj< Mat > K
 
SmartPetscObj< EPS > ePS
 
std::array< SmartPetscObj< Vec >, 6 > rigidBodyMotion
 
boost::shared_ptr< FEMethoddomianLhsFEPtr
 
boost::shared_ptr< FEMethoddomianRhsFEPtr
 
int fieldOrder = 2
 
SmartPetscObj< KSP > kspElastic
 
SmartPetscObj< DM > adjointDM
 
boost::shared_ptr< ObjectiveFunctionDatapythonPtr
 
std::vector< SmartPetscObj< Vec > > modeVecs
 
std::vector< std::array< double, 3 > > modeCentroids
 
std::vector< std::array< double, 6 > > modeBBoxes
 
SmartPetscObj< Vec > initialGeometry
 

Static Private Attributes

static std::vector< doublerZ
 
static std::vector< MatrixIntiI
 
static std::array< double, LAST_BBaveMaxMin
 
static int focalIndex
 
static int savitzkyGolayNormalisation
 
static const int * savitzkyGolayWeights
 
static ApproxFieldFunction< FIELD_DIMapproxFunction
 

Friends

struct TSPrePostProc
 

Detailed Description

Member Enumeration Documentation

◆ anonymous enum

anonymous enum
Enumerator
VOL 
COUNT 

Definition at line 222 of file plastic.cpp.

222{ VOL, COUNT };

◆ BoundingBox

enum Example::BoundingBox
private
Enumerator
CENTER_X 
CENTER_Y 
MAX_X 
MAX_Y 
MIN_X 
MIN_Y 
LAST_BB 
Examples
phase.cpp.

Definition at line 98 of file phase.cpp.

98 {
99 CENTER_X = 0,
100 CENTER_Y,
101 MAX_X,
102 MAX_Y,
103 MIN_X,
104 MIN_Y,
105 LAST_BB
106 };
@ MAX_X
Definition phase.cpp:101
@ MIN_X
Definition phase.cpp:103
@ MIN_Y
Definition phase.cpp:104
@ CENTER_X
Definition phase.cpp:99
@ MAX_Y
Definition phase.cpp:102
@ CENTER_Y
Definition phase.cpp:100
@ LAST_BB
Definition phase.cpp:105

Constructor & Destructor Documentation

◆ Example() [1/14]

Example::Example ( MoFEM::Interface & m_field)
inline

◆ Example() [2/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 420 of file dynamic_first_order_con_law.cpp.

420: mField(m_field) {}

◆ Example() [3/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 35 of file helmholtz.cpp.

35: mField(m_field) {}

◆ Example() [4/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 25 of file integration.cpp.

25: mField(m_field) {}

◆ Example() [5/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 50 of file plot_base.cpp.

50: mField(m_field) {}

◆ Example() [6/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 81 of file phase.cpp.

81: mField(m_field) {}

◆ Example() [7/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 53 of file approximaton.cpp.

53: mField(m_field) {}

◆ Example() [8/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 45 of file radiation.cpp.

45: mField(m_field) {}

◆ Example() [9/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 46 of file heat_method.cpp.

46: mField(m_field) {}

◆ Example() [10/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 80 of file elastic.cpp.

80: mField(m_field) {}

◆ Example() [11/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 49 of file eigen_elastic.cpp.

49: mField(m_field) {}

◆ Example() [12/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 52 of file nonlinear_dynamic_elastic.cpp.

52: mField(m_field) {}

◆ Example() [13/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 339 of file shallow_wave.cpp.

339: mField(m_field) {}

◆ Example() [14/14]

Example::Example ( MoFEM::Interface & m_field)
inline

Definition at line 130 of file adjoint.cpp.

130: mField(m_field) {}

Member Function Documentation

◆ assembleSystem() [1/10]

MoFEMErrorCode Example::assembleSystem ( )
private

[Push operators to pipeline]

[Adjoint modes]

[Boundary condition]

[Applying essential BC]

[Push operators to pipeline]

[Integration rule]

[Integration rule]

[Push domain stiffness matrix to pipeline]

[Push domain stiffness matrix to pipeline]

[Push Body forces]

[Push Body forces]

[Push natural boundary conditions]

[Push natural boundary conditions]

Examples
dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, mofem/tutorials/vec-7/adjoint.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 640 of file dynamic_first_order_con_law.cpp.

640 {
642 auto get_body_force = [this](const double, const double, const double) {
645 t_source(i) = 0.;
646 t_source(0) = 0.1;
647 t_source(1) = 1.;
648 return t_source;
649 };
650
651 // specific time scaling
652 auto get_time_scale = [this](const double time) {
653 return sin(time * omega * M_PI);
654 };
655
656 auto apply_rhs = [&](auto &pip) {
658
660 "GEOMETRY");
661
662 // Calculate Gradient of velocity
663 auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
665 "V", mat_v_grad_ptr));
666
667 auto gravity_vector_ptr = boost::make_shared<MatrixDouble>();
668 gravity_vector_ptr->resize(SPACE_DIM, 1);
669 auto set_body_force = [&]() {
672 auto t_force = getFTensor1FromMat<SPACE_DIM, 0>(*gravity_vector_ptr);
673 double unit_weight = 0.;
674 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-unit_weight", &unit_weight,
675 PETSC_NULLPTR);
676 t_force(i) = 0;
677 if (SPACE_DIM == 2) {
678 t_force(1) = -unit_weight;
679 } else if (SPACE_DIM == 3) {
680 t_force(2) = unit_weight;
681 }
683 };
684
685 CHKERR set_body_force();
686 pip.push_back(new OpBodyForce("V", gravity_vector_ptr,
687 [](double, double, double) { return 1.; }));
688
689 // Calculate unknown F
690 auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
692 "F", mat_H_tensor_ptr));
693
694 // // Calculate F
695 double tau = 0.2;
696 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-tau", &tau, PETSC_NULLPTR);
697
698 double xi = 0.;
699 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-xi", &xi, PETSC_NULLPTR);
700
701 // Calculate P stab
702 auto one = [&](const double, const double, const double) {
703 return 3. * bulk_modulus_K;
704 };
705 auto minus_one = [](const double, const double, const double) {
706 return -1.;
707 };
708
709 auto mat_dot_F_tensor_ptr = boost::make_shared<MatrixDouble>();
711 "F_dot", mat_dot_F_tensor_ptr));
712
713 // Calculate Gradient of Spatial Positions
714 auto mat_x_grad_ptr = boost::make_shared<MatrixDouble>();
716 "x_2", mat_x_grad_ptr));
717
718 auto mat_F_tensor_ptr = boost::make_shared<MatrixDouble>();
720 mat_F_tensor_ptr, mat_H_tensor_ptr));
721
722 auto mat_F_stab_ptr = boost::make_shared<MatrixDouble>();
724 mat_F_tensor_ptr, mat_F_stab_ptr, mat_dot_F_tensor_ptr, tau, xi,
725 mat_x_grad_ptr, mat_v_grad_ptr));
726
727 PetscBool is_linear_elasticity = PETSC_TRUE;
728 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_linear_elasticity",
729 &is_linear_elasticity, PETSC_NULLPTR);
730
731 auto mat_P_stab_ptr = boost::make_shared<MatrixDouble>();
732 if (is_linear_elasticity) {
735 mat_F_stab_ptr));
736 } else {
737 auto inv_F = boost::make_shared<MatrixDouble>();
738 auto det_ptr = boost::make_shared<VectorDouble>();
739
740 pip.push_back(
741 new OpInvertMatrix<SPACE_DIM>(mat_F_stab_ptr, det_ptr, inv_F));
742
743 // OpCalculatePiolaIncompressibleNH
746 mat_F_stab_ptr, inv_F, det_ptr));
747 }
748
749 pip.push_back(new OpGradTimesTensor2("V", mat_P_stab_ptr, minus_one));
750 pip.push_back(new OpRhsTestPiola("F", mat_v_grad_ptr, one));
751
753 };
754
755 auto *pipeline_mng = mField.getInterface<PipelineManager>();
756 CHKERR apply_rhs(pipeline_mng->getOpDomainExplicitRhsPipeline());
757
758 auto integration_rule = [](int, int, int approx_order) {
759 return 2 * approx_order;
760 };
761 CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);
762
764}
constexpr int SPACE_DIM
[Define dimension]
Definition adjoint.cpp:22
constexpr double shear_modulus_G
Definition adjoint.cpp:76
constexpr double bulk_modulus_K
Definition adjoint.cpp:75
@ H1
continuous field
Definition definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpBaseTimesVector< 1, SPACE_DIM *SPACE_DIM, 1 > OpRhsTestPiola
FormsIntegrators< DomainEleOp >::Assembly< AssemblyType::PETSC >::LinearForm< IntegrationType::GAUSS >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpGradTimesTensor2
constexpr double omega
Save field DOFS on vertices/tags.
auto integration_rule
FTensor::Index< 'i', SPACE_DIM > i
PetscErrorCode PetscOptionsGetReal(PetscOptions *, const char pre[], const char name[], PetscReal *dval, PetscBool *set)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
FTensor::Tensor1< FTensor::PackPtr< T *, S >, Tensor_Dim > getFTensor1FromMat(ublas::matrix< T, L, A > &data)
Get tensor rank 1 (vector) form data matrix.
static constexpr int approx_order
Add operators pushing bases from local to physical configuration.
Get values at integration pts for tensor field rank 2, i.e. matrix field.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
PipelineManager interface.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
DomainNaturalBCRhs::OpFlux< NaturalMeshsetType< BLOCKSET >, 1, SPACE_DIM > OpBodyForce

◆ assembleSystem() [2/10]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [3/10]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [4/10]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [5/10]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [6/10]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [7/10]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [8/10]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [9/10]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [10/10]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystemFlux()

MoFEMErrorCode Example::assembleSystemFlux ( )
private
Examples
phase.cpp.

◆ assembleSystemIntensity()

MoFEMErrorCode Example::assembleSystemIntensity ( )
private

[Calculate flux on boundary]

[Push operators to pipeline]

Examples
phase.cpp.

Definition at line 433 of file phase.cpp.

433 {
435
436 auto *pipeline_mng = mField.getInterface<PipelineManager>();
437
438 pipeline_mng->getDomainLhsFE().reset();
439 pipeline_mng->getDomainRhsFE().reset();
440 pipeline_mng->getBoundaryRhsFE().reset();
441
442 auto rule_vol = [](int, int, int order) { return 2 * (order + 1); };
443 pipeline_mng->setDomainLhsIntegrationRule(rule_vol);
444 pipeline_mng->setDomainRhsIntegrationRule(rule_vol);
445
446 auto det_ptr = boost::make_shared<VectorDouble>();
447 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
448 auto jac_ptr = boost::make_shared<MatrixDouble>();
449 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
450 pipeline_mng->getOpDomainLhsPipeline().push_back(
451 new OpCalculateHOJacForFace(jac_ptr));
452 pipeline_mng->getOpDomainLhsPipeline().push_back(
453 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
454 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMakeHdivFromHcurl());
455 pipeline_mng->getOpDomainLhsPipeline().push_back(
457 pipeline_mng->getOpDomainLhsPipeline().push_back(
458 new OpSetInvJacHcurlFace(inv_jac_ptr));
459 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
460
461 pipeline_mng->getOpDomainLhsPipeline().push_back(
462 new OpHdivHdiv("S", "S", lhsFlux));
463 auto unity = []() { return 1; };
464 pipeline_mng->getOpDomainLhsPipeline().push_back(
465 new OpHdivU("S", "PHI", unity, true));
466
467 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetHOWeightsOnFace());
468 pipeline_mng->getOpDomainRhsPipeline().push_back(
469 new OpDomainSource("PHI", rhsSource));
470
472}
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
constexpr int order
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)
Definition seepage.cpp:85
static double rhsSource(const double x, const double y, const double)
Definition phase.cpp:150
static double lhsFlux(const double x, const double y, const double)
Definition phase.cpp:165
Make Hdiv space from Hcurl space in 2d.
Modify integration weights on face to take in account higher-order geometry.
boost::shared_ptr< FEMethod > & getDomainLhsFE()

◆ bC() [1/2]

MoFEMErrorCode Example::bC ( )
private

[Create common data]

[Boundary condition]

Examples
plastic.cpp.

Definition at line 606 of file plastic.cpp.

606 {
608
610 auto bc_mng = mField.getInterface<BcManager>();
611
612 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
613 "U", 0, 0);
614 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
615 "U", 1, 1);
616 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
617 "U", 2, 2);
618 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
619 "REMOVE_ALL", "U", 0, 3);
620
621#ifdef ADD_CONTACT
622 for (auto b : {"FIX_X", "REMOVE_X"})
623 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
624 "SIGMA", 0, 0, false, true);
625 for (auto b : {"FIX_Y", "REMOVE_Y"})
626 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
627 "SIGMA", 1, 1, false, true);
628 for (auto b : {"FIX_Z", "REMOVE_Z"})
629 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
630 "SIGMA", 2, 2, false, true);
631 for (auto b : {"FIX_ALL", "REMOVE_ALL"})
632 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
633 "SIGMA", 0, 3, false, true);
634 CHKERR bc_mng->removeBlockDOFsOnEntities(
635 simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
636#endif
637
638 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
639 simple->getProblemName(), "U");
640
641 auto &bc_map = bc_mng->getBcMapByBlockName();
642 for (auto bc : bc_map)
643 MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
644
646}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69
#define MOFEM_LOG(channel, severity)
Log.
Simple interface for fast problem set-up.
Definition BcManager.hpp:29
Definition of the displacement bc data structure.
Definition BCData.hpp:72
Simple interface for fast problem set-up.
Definition Simple.hpp:27

◆ bC() [2/2]

MoFEMErrorCode Example::bC ( )
private

◆ boundaryCondition() [1/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

[Set up problem]

[Create common data]

[Boundary condition]

[Applying essential BC]

Examples
dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, mofem/tutorials/vec-7/adjoint.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 535 of file dynamic_first_order_con_law.cpp.

535 {
537
539 auto bc_mng = mField.getInterface<BcManager>();
540 auto *pipeline_mng = mField.getInterface<PipelineManager>();
541 auto time_scale = boost::make_shared<TimeScale>();
542
543 PetscBool sin_time_function = PETSC_FALSE;
544 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-sin_time_function",
545 &sin_time_function, PETSC_NULLPTR);
546
547 if (sin_time_function)
548 time_scale = boost::make_shared<DynamicFirstOrderConsSinusTimeScale>();
549 else
550 time_scale = boost::make_shared<DynamicFirstOrderConsConstantTimeScale>();
551
552 pipeline_mng->getBoundaryExplicitRhsFE().reset();
554 pipeline_mng->getOpBoundaryExplicitRhsPipeline(), {NOSPACE}, "GEOMETRY");
555
557 pipeline_mng->getOpBoundaryExplicitRhsPipeline(), mField, "V",
558 {time_scale}, "FORCE", Sev::inform);
559
560 auto integration_rule = [](int, int, int approx_order) {
561 return 2 * approx_order;
562 };
563
564 CHKERR pipeline_mng->setBoundaryExplicitRhsIntegrationRule(integration_rule);
565 CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);
566
567 CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
568 simple->getProblemName(), "V");
569
570 auto get_pre_proc_hook = [&]() {
572 mField, pipeline_mng->getDomainExplicitRhsFE(), {time_scale});
573 };
574 pipeline_mng->getDomainExplicitRhsFE()->preProcessHook = get_pre_proc_hook();
575
577}
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25

◆ boundaryCondition() [2/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [3/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [4/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [5/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [6/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [7/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [8/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [9/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [10/10]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ calculateFlux()

MoFEMErrorCode Example::calculateFlux ( double & calc_flux)
private

[Set up problem]

[Calculate flux on boundary]

Examples
phase.cpp.

Definition at line 402 of file phase.cpp.

402 {
404 auto pipeline_mng = mField.getInterface<PipelineManager>();
405
406 pipeline_mng->getDomainLhsFE().reset();
407 pipeline_mng->getDomainRhsFE().reset();
408 pipeline_mng->getBoundaryRhsFE().reset();
409
410 auto rule_vol = [](int, int, int order) { return 2 * (order + 1); };
411 pipeline_mng->setBoundaryRhsIntegrationRule(rule_vol);
412
413 auto flux_ptr = boost::make_shared<MatrixDouble>();
414 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
416 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
417 new OpCalculateHVecVectorField<3>("S", flux_ptr));
418 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
419 new BoundaryOp(flux_ptr, calc_flux));
420
421 calc_flux = 0;
422 CHKERR pipeline_mng->loopFiniteElements();
423 double global_flux_assembeld = 0;
424 MPI_Allreduce(&calc_flux, &global_flux_assembeld, 1, MPI_DOUBLE, MPI_SUM,
425 mField.get_comm());
426 calc_flux = global_flux_assembeld;
427
429}
virtual MPI_Comm & get_comm() const =0
Get vector field for H-div approximation.

◆ calculateGradient()

MoFEMErrorCode Example::calculateGradient ( PetscReal * objective_function_value,
Vec objective_function_gradient,
Vec adjoint_vector )
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 1716 of file adjoint.cpp.

1718 {
1719 MOFEM_LOG_CHANNEL("WORLD");
1721 auto simple = mField.getInterface<Simple>();
1722
1723 auto ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
1724
1725 auto get_essential_fe = [this]() {
1726 auto post_proc_rhs = boost::make_shared<FEMethod>();
1727 auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
1729
1731 post_proc_rhs, 0)();
1733 };
1734 post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
1735 return post_proc_rhs;
1736 };
1737
1738 auto get_fd_direvative_fe = [&]() {
1739 auto fe = boost::make_shared<DomainEle>(mField);
1740 fe->getRuleHook = [](int, int, int p_data) {
1741 return 2 * p_data + p_data - 1;
1742 };
1743 auto &pip = fe->getOpPtrVector();
1745 // Add RHS operators for internal forces
1746 constexpr bool debug = false;
1747 if constexpr (debug) {
1749 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
1750 pip.push_back(
1751 new OpAdJointTestOp<SPACE_DIM, I, DomainBaseOp>("U", common_ptr));
1752 } else {
1754 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
1755 }
1756 return fe;
1757 };
1758
1759 auto calulate_fd_residual = [&](auto eps, auto diff_vec, auto fd_vec) {
1761
1762 constexpr bool debug = false;
1763
1764 auto geom_norm = [](MoFEM::Interface &mField) {
1766 auto field_blas = mField.getInterface<FieldBlas>();
1767 double nrm2 = 0.0;
1768 auto norm2_field = [&](const double val) {
1769 nrm2 += val * val;
1770 return val;
1771 };
1772 CHKERR field_blas->fieldLambdaOnValues(norm2_field, "GEOMETRY");
1773 MPI_Allreduce(MPI_IN_PLACE, &nrm2, 1, MPI_DOUBLE, MPI_SUM,
1774 mField.get_comm());
1775 MOFEM_LOG("WORLD", Sev::inform) << "Geometry norm: " << sqrt(nrm2);
1777 };
1778
1779 if constexpr (debug)
1780 CHKERR geom_norm(mField);
1781
1782 auto initial_current_geometry = createDMVector(adjointDM);
1783 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1784 "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
1785 initial_current_geometry, INSERT_VALUES, SCATTER_FORWARD);
1786 CHKERR VecAssemblyBegin(initial_current_geometry);
1787 CHKERR VecAssemblyEnd(initial_current_geometry);
1788
1789 if constexpr (debug)
1790 CHKERR geom_norm(mField);
1791
1792 auto perturb_geometry = [&](auto eps, auto diff_vec) {
1794 auto current_geometry = vectorDuplicate(initial_current_geometry);
1795 CHKERR VecCopy(initial_current_geometry, current_geometry);
1796 CHKERR VecAXPY(current_geometry, eps, diff_vec);
1797 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1798 "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
1799 current_geometry, INSERT_VALUES, SCATTER_REVERSE);
1801 };
1802
1803 auto fe = get_fd_direvative_fe();
1804 auto fp = vectorDuplicate(diff_vec);
1805 auto fm = vectorDuplicate(diff_vec);
1806 auto calc_impl = [&](auto f, auto eps) {
1808 CHKERR VecZeroEntries(f);
1809 fe->f = f;
1810 CHKERR perturb_geometry(eps, diff_vec);
1812 simple->getDomainFEName(), fe);
1813 CHKERR VecAssemblyBegin(f);
1814 CHKERR VecAssemblyEnd(f);
1815 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
1816 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
1817 auto post_proc_rhs = get_essential_fe();
1818 post_proc_rhs->f = f;
1820 post_proc_rhs.get());
1822 };
1823 CHKERR calc_impl(fp, eps);
1824 CHKERR calc_impl(fm, -eps);
1825 CHKERR VecWAXPY(fd_vec, -1.0, fm, fp);
1826 CHKERR VecScale(fd_vec, 1.0 / (2.0 * eps));
1827
1828 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1829 "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
1830 initial_current_geometry, INSERT_VALUES, SCATTER_REVERSE);
1831
1832 if constexpr (debug)
1833 CHKERR geom_norm(mField);
1834
1836 };
1837
1838 auto get_direvative_fe = [&](auto diff_vec) {
1839 auto fe_adjoint = boost::make_shared<DomainEle>(mField);
1840 fe_adjoint->getRuleHook = [](int, int, int p_data) {
1841 return 2 * p_data + p_data - 1;
1842 };
1843 auto &pip = fe_adjoint->getOpPtrVector();
1844
1845 auto jac_ptr = boost::make_shared<MatrixDouble>();
1846 auto det_ptr = boost::make_shared<VectorDouble>();
1847 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1848 auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
1849 auto cof_ptr = boost::make_shared<VectorDouble>();
1850
1851 using OpCoFactor =
1852 AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
1853
1855
1857 "GEOMETRY", jac_ptr));
1859 "U", diff_jac_ptr,
1860 diff_vec)); // Note: that vector is stored on displacemnt vector, that
1861 // why is used here
1862
1863 pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
1864
1865 // Add RHS operators for internal forces
1867 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
1869 "U", common_ptr, jac_ptr, diff_jac_ptr, cof_ptr));
1870
1871 return fe_adjoint;
1872 };
1873
1874 auto get_objective_fe = [&](auto diff_vec, auto grad_vec,
1875 auto glob_objective_ptr,
1876 auto glob_objective_grad_ptr) {
1877 auto fe_adjoint = boost::make_shared<DomainEle>(mField);
1878 fe_adjoint->getRuleHook = [](int, int, int p_data) {
1879 return 2 * p_data + p_data - 1;
1880 };
1881 auto &pip = fe_adjoint->getOpPtrVector();
1883
1884 auto jac_ptr = boost::make_shared<MatrixDouble>();
1885 auto det_ptr = boost::make_shared<VectorDouble>();
1886 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1887 auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
1888 auto cof_ptr = boost::make_shared<VectorDouble>();
1889 auto d_grad_ptr = boost::make_shared<MatrixDouble>();
1890 auto u_ptr = boost::make_shared<MatrixDouble>();
1891
1892 using OpCoFactor =
1893 AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
1895 "GEOMETRY", jac_ptr));
1897 "U", diff_jac_ptr,
1898 diff_vec)); // Note: that vector is stored on displacemnt vector, that
1899 // why is used here
1900 pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
1902 "U", d_grad_ptr, grad_vec));
1903 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1904
1906 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
1907 pip.push_back(new OpAdJointObjective(
1908 pythonPtr, common_ptr, jac_ptr, diff_jac_ptr, cof_ptr, d_grad_ptr,
1909 u_ptr, glob_objective_ptr, glob_objective_grad_ptr));
1910
1911 return fe_adjoint;
1912 };
1913
1914 auto dm = simple->getDM();
1915 auto f = createDMVector(dm);
1916 auto d = vectorDuplicate(f);
1917 auto dm_diff_vec = vectorDuplicate(d);
1918
1919 auto adjoint_fe = get_direvative_fe(dm_diff_vec);
1920 auto glob_objective_ptr = boost::make_shared<double>(0.0);
1921 auto glob_objective_grad_ptr = boost::make_shared<double>(0.0);
1922 auto objective_fe = get_objective_fe(dm_diff_vec, d, glob_objective_ptr,
1923 glob_objective_grad_ptr);
1924
1925 auto set_variance_of_geometry = [&](auto mode, auto mod_vec) {
1927 CHKERR DMoFEMMeshToLocalVector(adjointDM, mod_vec, INSERT_VALUES,
1928 SCATTER_REVERSE);
1929 CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
1930 simple->getProblemName(), "U", "ADJOINT_FIELD", RowColData::ROW,
1931 dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
1932 CHKERR VecGhostUpdateBegin(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
1933 CHKERR VecGhostUpdateEnd(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
1935 };
1936
1937 auto calculate_variance_internal_forces = [&](auto mode, auto mod_vec) {
1939 CHKERR VecZeroEntries(f);
1940 CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
1941 CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
1942 adjoint_fe->f = f;
1943 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), adjoint_fe);
1944 CHKERR VecAssemblyBegin(f);
1945 CHKERR VecAssemblyEnd(f);
1946 CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
1947 CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
1948 auto post_proc_rhs = get_essential_fe();
1949 post_proc_rhs->f = f;
1951 post_proc_rhs.get());
1952 CHKERR VecScale(f, -1.0);
1953
1954#ifndef NDEBUG
1955 constexpr bool debug = false;
1956 if constexpr (debug) {
1957 double norm0;
1958 CHKERR VecNorm(f, NORM_2, &norm0);
1959 auto fd_check = vectorDuplicate(f);
1960 double eps = 1e-5;
1961 CHKERR calulate_fd_residual(eps, dm_diff_vec, fd_check);
1962 double nrm;
1963 CHKERR VecAXPY(fd_check, -1.0, f);
1964 CHKERR VecNorm(fd_check, NORM_2, &nrm);
1965 MOFEM_LOG("WORLD", Sev::inform)
1966 << " FD check for internal forces [ " << mode << " ]: " << nrm
1967 << " / " << norm0 << " ( " << (nrm / norm0) << " )";
1968 }
1969#endif
1970
1972 };
1973
1974 auto calulate_variance_of_displacemnte = [&](auto mode, auto mod_vec) {
1976 CHKERR KSPSolve(kspElastic, f, d);
1977 CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
1978 CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
1980 };
1981
1982 auto calculate_variance_of_objective_function = [&](auto mode, auto mod_vec) {
1984 *glob_objective_ptr = 0.0;
1985 *glob_objective_grad_ptr = 0.0;
1986 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
1987 objective_fe);
1988 CHKERR VecSetValue(objective_function_gradient, mode,
1989 *glob_objective_grad_ptr, ADD_VALUES);
1990 std::array<double, 2> array = {*glob_objective_ptr,
1991 *glob_objective_grad_ptr};
1992 MPI_Allreduce(MPI_IN_PLACE, array.data(), 2, MPI_DOUBLE, MPI_SUM,
1993 mField.get_comm());
1994 *glob_objective_ptr = array[0];
1995 *glob_objective_grad_ptr = array[1];
1996 (*objective_function_value) += *glob_objective_ptr;
1997 MOFEM_LOG_CHANNEL("WORLD");
1998 MOFEM_LOG("WORLD", Sev::verbose) << " Objective gradient [ " << mode
1999 << " ]: " << *glob_objective_grad_ptr;
2001 };
2002
2003 CHKERR VecZeroEntries(objective_function_gradient);
2004 CHKERR VecZeroEntries(adjoint_vector);
2005 *objective_function_value = 0.0;
2006 int mode = 0;
2007 for (auto mod_vec : modeVecs) {
2008
2009 CHKERR set_variance_of_geometry(mode, mod_vec);
2010 CHKERR calculate_variance_internal_forces(mode, mod_vec);
2011 CHKERR calulate_variance_of_displacemnte(mode, mod_vec);
2012 CHKERR calculate_variance_of_objective_function(mode, mod_vec);
2013
2014 CHKERR VecAXPY(adjoint_vector, *glob_objective_grad_ptr, dm_diff_vec);
2015 ++mode;
2016 }
2017
2018 (*objective_function_value) /= modeVecs.size();
2019 MOFEM_LOG("WORLD", Sev::verbose)
2020 << "Objective function: " << *glob_objective_ptr;
2021
2022 CHKERR VecAssemblyBegin(objective_function_gradient);
2023 CHKERR VecAssemblyEnd(objective_function_gradient);
2024
2025 CHKERR VecAssemblyBegin(adjoint_vector);
2026 CHKERR VecAssemblyEnd(adjoint_vector);
2027 CHKERR VecGhostUpdateBegin(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
2028 CHKERR VecGhostUpdateEnd(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
2029
2031}
Range get_range_from_block(MoFEM::Interface &m_field, const std::string block_name, int dim)
Definition adjoint.cpp:2579
static const double eps
@ ROW
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
auto createDMVector(DM dm)
Get smart vector from DM.
Definition DMMoFEM.hpp:1102
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
auto commonDataFactory(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, std::string block_name, Sev sev, double scale=1)
Definition HookeOps.hpp:208
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev, bool is_non_linear=false)
Factory function to create and push internal force operators for the domain RHS. (RHS first overload)...
Definition HookeOps.hpp:242
static const bool debug
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
SmartPetscObj< KSP > kspElastic
Definition adjoint.cpp:155
std::vector< SmartPetscObj< Vec > > modeVecs
Definition adjoint.cpp:158
boost::shared_ptr< ObjectiveFunctionData > pythonPtr
Definition adjoint.cpp:157
SmartPetscObj< DM > adjointDM
Definition adjoint.cpp:156
Deprecated interface functions.
Class (Function) to enforce essential constrains on the right hand side diagonal.
Definition Essential.hpp:41
Basic algebra on fields.
Definition FieldBlas.hpp:21
MoFEMErrorCode fieldLambdaOnValues(OneFieldFunctionOnValues lambda, const std::string field_name, Range *ents_ptr=nullptr)
field lambda
Definition FieldBlas.cpp:21
Get values at integration pts for tensor field rank 1, i.e. vector field.
Vector manager is used to create vectors \mofem_vectors.

◆ checkResults() [1/11]

MoFEMErrorCode Example::checkResults ( )
private

[Postprocess results]

[Print results]

[Check]

[Check results]

[Test example]

Examples
dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 1182 of file dynamic_first_order_con_law.cpp.

1182 {
1185}

◆ checkResults() [2/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ checkResults() [3/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ checkResults() [4/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ checkResults() [5/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ checkResults() [6/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ checkResults() [7/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ checkResults() [8/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ checkResults() [9/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ checkResults() [10/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ checkResults() [11/11]

MoFEMErrorCode Example::checkResults ( )
private

◆ createCommonData() [1/8]

MoFEMErrorCode Example::createCommonData ( )
private

[Set up problem]

[Set integration rule]

[Create common data]

< true if tau order is set

< true if tau order is set

Examples
eigen_elastic.cpp, heat_method.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 480 of file plastic.cpp.

480 {
482
483 auto get_command_line_parameters = [&]() {
485
486 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-scale", &scale, PETSC_NULLPTR);
487 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
488 &young_modulus, PETSC_NULLPTR);
489 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
490 &poisson_ratio, PETSC_NULLPTR);
491 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-hardening", &H, PETSC_NULLPTR);
492 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-hardening_viscous", &visH,
493 PETSC_NULLPTR);
494 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-yield_stress", &sigmaY,
495 PETSC_NULLPTR);
496 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn0", &cn0, PETSC_NULLPTR);
497 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn1", &cn1, PETSC_NULLPTR);
498 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-zeta", &zeta, PETSC_NULLPTR);
499 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-Qinf", &Qinf, PETSC_NULLPTR);
500 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-b_iso", &b_iso, PETSC_NULLPTR);
501 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-C1_k", &C1_k, PETSC_NULLPTR);
502 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-large_strains",
503 &is_large_strains, PETSC_NULLPTR);
504 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-set_timer", &set_timer,
505 PETSC_NULLPTR);
506 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
507 PETSC_NULLPTR);
508
509 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
510 PetscBool tau_order_is_set; ///< true if tau order is set
511 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-tau_order", &tau_order,
512 &tau_order_is_set);
513 PetscBool ep_order_is_set; ///< true if tau order is set
514 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-ep_order", &ep_order,
515 &ep_order_is_set);
516 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
517 PETSC_NULLPTR);
518
519 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho", &rho, PETSC_NULLPTR);
520 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_damping",
521 &alpha_damping, PETSC_NULLPTR);
522
523 MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;
524 MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;
525 MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;
526 MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;
527 MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;
528 MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;
529 MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;
530 MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;
531 MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;
532 MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;
533 MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;
534
535 if (tau_order_is_set == PETSC_FALSE)
536 tau_order = order - 2;
537 if (ep_order_is_set == PETSC_FALSE)
538 ep_order = order - 1;
539
540 MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
541 MOFEM_LOG("PLASTICITY", Sev::inform)
542 << "Ep approximation order " << ep_order;
543 MOFEM_LOG("PLASTICITY", Sev::inform)
544 << "Tau approximation order " << tau_order;
545 MOFEM_LOG("PLASTICITY", Sev::inform)
546 << "Geometry approximation order " << geom_order;
547
548 MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;
549 MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;
550
551 PetscBool is_scale = PETSC_TRUE;
552 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_scale", &is_scale,
553 PETSC_NULLPTR);
554 if (is_scale) {
556 }
557
558 MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
559
560#ifdef ADD_CONTACT
561 CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn_contact",
562 &ContactOps::cn_contact, PETSC_NULLPTR);
563 MOFEM_LOG("CONTACT", Sev::inform)
564 << "cn_contact " << ContactOps::cn_contact;
565#endif // ADD_CONTACT
566
567 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-quasi_static",
568 &is_quasi_static, PETSC_NULLPTR);
569 MOFEM_LOG("PLASTICITY", Sev::inform)
570 << "Is quasi static: " << (is_quasi_static ? "true" : "false");
571
573 };
574
575 CHKERR get_command_line_parameters();
576
577#ifdef ADD_CONTACT
578 #ifdef ENABLE_PYTHON_BINDING
579 auto file_exists = [](std::string myfile) {
580 std::ifstream file(myfile.c_str());
581 if (file) {
582 return true;
583 }
584 return false;
585 };
586 char sdf_file_name[255] = "sdf.py";
587 CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-sdf_file",
588 sdf_file_name, 255, PETSC_NULLPTR);
589
590 if (file_exists(sdf_file_name)) {
591 MOFEM_LOG("CONTACT", Sev::inform) << sdf_file_name << " file found";
592 sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
593 CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
594 ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
595 } else {
596 MOFEM_LOG("CONTACT", Sev::warning) << sdf_file_name << " file NOT found";
597 }
598 #endif
599#endif // ADD_CONTACT
600
602}
constexpr double poisson_ratio
Definition adjoint.cpp:74
constexpr double young_modulus
Definition adjoint.cpp:73
double cn_contact
Definition contact.cpp:99
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
double C1_k
Kinematic hardening.
Definition plastic.cpp:133
double Qinf
Saturation yield stress.
Definition plastic.cpp:131
double rho
Definition plastic.cpp:144
int atom_test
Atom test.
Definition plastic.cpp:121
PetscBool is_quasi_static
Definition plastic.cpp:143
double alpha_damping
Definition plastic.cpp:145
double visH
Viscous hardening.
Definition plastic.cpp:129
PetscBool set_timer
Set timer.
Definition plastic.cpp:118
double scale
Definition plastic.cpp:123
double zeta
Viscous hardening.
Definition plastic.cpp:130
double H
Hardening.
Definition plastic.cpp:128
int tau_order
Order of tau files.
Definition plastic.cpp:139
double cn0
Definition plastic.cpp:135
double b_iso
Saturation exponent.
Definition plastic.cpp:132
PetscBool is_large_strains
Large strains.
Definition plastic.cpp:117
int geom_order
Order if fixed.
Definition plastic.cpp:141
double sigmaY
Yield stress.
Definition plastic.cpp:127
int ep_order
Order of ep files.
Definition plastic.cpp:140
double cn1
Definition plastic.cpp:136

◆ createCommonData() [2/8]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [3/8]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [4/8]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [5/8]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [6/8]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [7/8]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [8/8]

MoFEMErrorCode Example::createCommonData ( )
private

◆ getCoordsInImage()

std::pair< int, int > Example::getCoordsInImage ( double x,
double y )
staticprivate
Examples
phase.cpp.

Definition at line 131 of file phase.cpp.

131 {
132
133 auto &m = iI[focalIndex];
134 x -= aveMaxMin[MIN_X];
135 y -= aveMaxMin[MIN_Y];
136 x *= (m.size1() - 1) / (aveMaxMin[MAX_X] - aveMaxMin[MIN_X]);
137 y *= (m.size2() - 1) / (aveMaxMin[MAX_Y] - aveMaxMin[MIN_Y]);
138 const auto p = std::make_pair<int, int>(std::round(x), std::round(y));
139
140#ifndef NDEBUG
141 if (p.first < 0 && p.first >= m.size1())
142 THROW_MESSAGE("Wrong index");
143 if (p.second < 0 && p.second >= m.size2())
144 THROW_MESSAGE("Wrong index");
145#endif
146
147 return p;
148}
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
FTensor::Index< 'm', 3 > m
static std::array< double, LAST_BB > aveMaxMin
Definition phase.cpp:110
static int focalIndex
Definition phase.cpp:112
static std::vector< MatrixInt > iI
Definition phase.cpp:109

◆ integrateElements()

MoFEMErrorCode Example::integrateElements ( )
private

[Push operators to pipeline]

[Integrate]

Definition at line 218 of file integration.cpp.

218 {
220 // Zero global vector
221 CHKERR VecZeroEntries(commonDataPtr->petscVec);
222
223 // Integrate elements by executing operators in the pipeline
225 CHKERR pipeline_mng->loopFiniteElements();
226
227 // Assemble MPI vector
228 CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
229 CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
231}
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::shared_ptr< CommonData > commonDataPtr

◆ integrationRule()

static int Example::integrationRule ( int ,
int ,
int p_data )
inlinestaticprivate

Definition at line 52 of file radiation.cpp.

52{ return 2 * p_data; };

◆ kspSolve()

MoFEMErrorCode Example::kspSolve ( )
private

[Push operators to pipeline]

[Solve]

Definition at line 230 of file radiation.cpp.

230 {
234 auto ts = pipeline_mng->createTSIM();
235
236 double ftime = 1;
237 CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
238 CHKERR TSSetFromOptions(ts);
239 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
240
241 auto T = createDMVector(simple->getDM());
242 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
243 SCATTER_FORWARD);
244
245 CHKERR TSSolve(ts, T);
246 CHKERR TSGetTime(ts, &ftime);
247
248 PetscInt steps, snesfails, rejects, nonlinits, linits;
249 CHKERR TSGetTimeStepNumber(ts, &steps);
250 CHKERR TSGetSNESFailures(ts, &snesfails);
251 CHKERR TSGetStepRejections(ts, &rejects);
252 CHKERR TSGetSNESIterations(ts, &nonlinits);
253 CHKERR TSGetKSPIterations(ts, &linits);
254 MOFEM_LOG_C("EXAMPLE", Sev::inform,
255 "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
256 "%d, linits %d",
257 steps, rejects, snesfails, ftime, nonlinits, linits);
258
260}
#define MOFEM_LOG_C(channel, severity, format,...)
SmartPetscObj< TS > createTSIM(SmartPetscObj< DM > dm=nullptr)
Create TS (time) implicit solver.

◆ lhsFlux()

double Example::lhsFlux ( const double x,
const double y,
const double  )
staticprivate
Examples
phase.cpp.

Definition at line 165 of file phase.cpp.

165 {
166 const auto idx = getCoordsInImage(x, y);
167 const auto &m = iI[focalIndex];
168 return 1. / m(idx.first, idx.second);
169}
static std::pair< int, int > getCoordsInImage(double x, double y)
Definition phase.cpp:131

◆ OPs() [1/2]

MoFEMErrorCode Example::OPs ( )
private

[Boundary condition]

[Push operators to pipeline]

[Only used for dynamics]

[Only used for dynamics]

[Only used for dynamics]

[Only used for dynamics]

Examples
plastic.cpp.

Definition at line 650 of file plastic.cpp.

650 {
652 auto pip_mng = mField.getInterface<PipelineManager>();
653
654 auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
655
656 auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
657
658 auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {
660
662 pip, {HDIV}, "GEOMETRY");
663 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
664
665 // Add Natural BCs to LHS
667 pip, mField, "U", Sev::inform);
668
669#ifdef ADD_CONTACT
671 CHKERR
673 pip, "SIGMA", "U");
674 CHKERR
676 mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
677 vol_rule);
678#endif // ADD_CONTACT
679
681 };
682
683 auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
685
687 pip, {HDIV}, "GEOMETRY");
688 pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
689
690 // Add Natural BCs to RHS
692 pip, mField, "U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
693
694#ifdef ADD_CONTACT
696 pip, "SIGMA", "U");
697#endif // ADD_CONTACT
698
700 };
701
702 auto add_domain_ops_lhs = [this](auto &pip) {
705 pip, {H1, HDIV}, "GEOMETRY");
706
707 if (is_quasi_static == PETSC_FALSE) {
708
709 //! [Only used for dynamics]
712 //! [Only used for dynamics]
713
714 auto get_inertia_and_mass_damping = [this](const double, const double,
715 const double) {
716 auto *pip = mField.getInterface<PipelineManager>();
717 auto &fe_domain_lhs = pip->getDomainLhsFE();
718 return (rho / scale) * fe_domain_lhs->ts_aa +
719 (alpha_damping / scale) * fe_domain_lhs->ts_a;
720 };
721 pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
722 }
723
725 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
726
728 };
729
730 auto add_domain_ops_rhs = [this](auto &pip) {
732
734 pip, {H1, HDIV}, "GEOMETRY");
735
737 pip, mField, "U",
738 {boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
739 Sev::inform);
740
741 // only in case of dynamics
742 if (is_quasi_static == PETSC_FALSE) {
743
744 //! [Only used for dynamics]
747 //! [Only used for dynamics]
748
749 auto mat_acceleration = boost::make_shared<MatrixDouble>();
751 "U", mat_acceleration));
752 pip.push_back(
753 new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
754 return rho / scale;
755 }));
756 if (alpha_damping > 0) {
757 auto mat_velocity = boost::make_shared<MatrixDouble>();
758 pip.push_back(
759 new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
760 pip.push_back(
761 new OpInertiaForce("U", mat_velocity, [](double, double, double) {
762 return alpha_damping / scale;
763 }));
764 }
765 }
766
768 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
769
770#ifdef ADD_CONTACT
772 pip, "SIGMA", "U");
773#endif // ADD_CONTACT
774
776 };
777
778 CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
779 CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
780
781 // Boundary
782 CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
783 CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
784
785 CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
786 CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
787
788 CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
789 CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
790
791 auto create_reaction_pipeline = [&](auto &pip) {
794 pip, {H1}, "GEOMETRY");
796 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
798 };
799
800 reactionFe = boost::make_shared<DomainEle>(mField);
801 reactionFe->getRuleHook = vol_rule;
802 CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
803 reactionFe->postProcessHook =
805
807}
@ HDIV
field with continuous normal traction
Definition definitions.h:87
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesVector< 1, SPACE_DIM, 1 > OpInertiaForce
MoFEMErrorCode opFactoryDomainRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
MoFEMErrorCode opFactoryBoundaryToDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string fe_domain_name, std::string sigma, std::string u, std::string geom, ForcesAndSourcesCore::RuleHookFun rule, bool is_axisymmetric=false)
MoFEMErrorCode opFactoryBoundaryRhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
MoFEMErrorCode opFactoryBoundaryLhs(boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string sigma, std::string u, bool is_axisymmetric=false)
MoFEMErrorCode opFactoryDomainReactions(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
MoFEMErrorCode opFactoryDomainRhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, std::string block_name, Pip &pip, std::string u, std::string ep, std::string tau)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition seepage.cpp:55
boost::shared_ptr< DomainEle > reactionFe
Definition plastic.cpp:235
Class (Function) to calculate residual side diagonal.
Definition Essential.hpp:49
Approximate field values for given petsc vector.
constexpr AssemblyType AT

◆ OPs() [2/2]

MoFEMErrorCode Example::OPs ( )
private

◆ outputResults() [1/10]

MoFEMErrorCode Example::outputResults ( )
private

[Solve]

[Postprocess results]

[Postprocess clean]

[Postprocess clean]

[Postprocess initialise]

[Postprocess initialise]

Examples
dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, phase.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 1160 of file dynamic_first_order_con_law.cpp.

1160 {
1162 PetscBool test_flg = PETSC_FALSE;
1163 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test", &test_flg, PETSC_NULLPTR);
1164 if (test_flg) {
1165 auto *simple = mField.getInterface<Simple>();
1166 auto T = createDMVector(simple->getDM());
1167 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
1168 SCATTER_FORWARD);
1169 double nrm2;
1170 CHKERR VecNorm(T, NORM_2, &nrm2);
1171 MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;
1172 constexpr double regression_value = 0.0194561;
1173 if (fabs(nrm2 - regression_value) > 1e-2)
1174 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1175 "Regression test failed; wrong norm value.");
1176 }
1178}
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40

◆ outputResults() [2/10]

MoFEMErrorCode Example::outputResults ( )
private

◆ outputResults() [3/10]

MoFEMErrorCode Example::outputResults ( )
private

◆ outputResults() [4/10]

MoFEMErrorCode Example::outputResults ( )
private

◆ outputResults() [5/10]

MoFEMErrorCode Example::outputResults ( )
private

◆ outputResults() [6/10]

MoFEMErrorCode Example::outputResults ( )
private

◆ outputResults() [7/10]

MoFEMErrorCode Example::outputResults ( )
private

◆ outputResults() [8/10]

MoFEMErrorCode Example::outputResults ( )
private

◆ outputResults() [9/10]

MoFEMErrorCode Example::outputResults ( )
private

◆ outputResults() [10/10]

MoFEMErrorCode Example::outputResults ( const int i)
private

[Solve]

[Postprocess results]

Definition at line 504 of file phase.cpp.

504 {
507 pipeline_mng->getDomainLhsFE().reset();
508 pipeline_mng->getDomainRhsFE().reset();
509 pipeline_mng->getBoundaryRhsFE().reset();
510
512
513 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
514 auto jac_ptr = boost::make_shared<MatrixDouble>();
515 post_proc_fe->getOpPtrVector().push_back(
516 new OpCalculateHOJacForFace(jac_ptr));
517 post_proc_fe->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
518 post_proc_fe->getOpPtrVector().push_back(
520
521 auto s_ptr = boost::make_shared<VectorDouble>();
522 auto phi_ptr = boost::make_shared<MatrixDouble>();
523 post_proc_fe->getOpPtrVector().push_back(
524 new OpCalculateScalarFieldValues("S", s_ptr));
525 post_proc_fe->getOpPtrVector().push_back(
526 new OpCalculateHVecVectorField<3>("PHI", phi_ptr));
527
529
530 post_proc_fe->getOpPtrVector().push_back(
531
532 new OpPPMap(post_proc_fe->getPostProcMesh(),
533 post_proc_fe->getMapGaussPts(),
534
535 OpPPMap::DataMapVec{{"S", s_ptr}},
536
537 OpPPMap::DataMapMat{{"PHI", phi_ptr}},
538
540
542
543 )
544
545 );
546
547 // post_proc_fe->addFieldValuesPostProc("S");
548 // post_proc_fe->addFieldValuesPostProc("PHI");
549
550
551 pipeline_mng->getDomainRhsFE() = post_proc_fe;
552 CHKERR pipeline_mng->loopFiniteElements();
553 CHKERR post_proc_fe->writeFile("out_" + boost::lexical_cast<std::string>(i) +
554 ".h5m");
556}
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Get value at integration points for scalar 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
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()

◆ postProcess() [1/2]

MoFEMErrorCode Example::postProcess ( )
private

[Integrate]

[Solve]

[Print results]

[Postprocess results]

Definition at line 235 of file integration.cpp.

235 {
237 const double *array;
238 CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
239 if (mField.get_comm_rank() == 0) {
240 MOFEM_LOG_C("SELF", Sev::inform, "Mass %6.4e", array[CommonData::ZERO]);
241 MOFEM_LOG_C("SELF", Sev::inform,
242 "First moment of inertia [ %6.4e, %6.4e, %6.4e ]",
244 array[CommonData::FIRST_Z]);
245 MOFEM_LOG_C("SELF", Sev::inform,
246 "Second moment of inertia [ %6.4e, %6.4e, %6.4e; %6.4e %6.4e; "
247 "%6.4e ]",
251 }
252 CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
254}
virtual int get_comm_rank() const =0

◆ postProcess() [2/2]

MoFEMErrorCode Example::postProcess ( )
private

◆ postprocessElastic()

MoFEMErrorCode Example::postprocessElastic ( int iter,
SmartPetscObj< Vec > adjoint_vector = nullptr )
private

[Solve]

[Postprocess results]

[Postprocess clean]

[Postprocess clean]

[Postprocess initialise]

[Postprocess initialise]

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 1164 of file adjoint.cpp.

1165 {
1167 auto simple = mField.getInterface<Simple>();
1168 auto pip = mField.getInterface<PipelineManager>();
1169 auto det_ptr = boost::make_shared<VectorDouble>();
1170 auto jac_ptr = boost::make_shared<MatrixDouble>();
1171 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
1172 //! [Postprocess clean]
1173 pip->getDomainRhsFE().reset();
1174 pip->getDomainLhsFE().reset();
1175 pip->getBoundaryRhsFE().reset();
1176 pip->getBoundaryLhsFE().reset();
1177 //! [Postprocess clean]
1178
1179 //! [Postprocess initialise]
1180 auto post_proc_mesh = boost::make_shared<moab::Core>();
1181 auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
1182 mField, post_proc_mesh);
1183 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
1184 mField, post_proc_mesh);
1185 //! [Postprocess initialise]
1186
1187 // lambada function to calculate displacement, strain and stress
1188 auto calculate_stress_ops = [&](auto &pip) {
1189 // Add H1 geometry ops
1191
1192 // Use HookeOps commonDataFactory to get matStrainPtr and matCauchyStressPtr
1194 mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
1195
1196 // Store U and GEOMETRY values if needed
1197 auto u_ptr = boost::make_shared<MatrixDouble>();
1198 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
1199 auto x_ptr = boost::make_shared<MatrixDouble>();
1200 pip.push_back(
1201 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
1202
1203 using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
1204
1205 DataMapMat vec_map = {{"U", u_ptr}, {"GEOMETRY", x_ptr}};
1206 DataMapMat sym_map = {{"STRAIN", common_ptr->getMatStrain()},
1207 {"STRESS", common_ptr->getMatCauchyStress()}};
1208
1209 if (adjoint_vector) {
1210 auto adjoint_ptr = boost::make_shared<MatrixDouble>();
1212 "U", adjoint_ptr, adjoint_vector));
1213 vec_map["ADJOINT"] = adjoint_ptr;
1214 }
1215
1216 // Return what you need: displacements, coordinates, strain, stress
1217 return boost::make_tuple(u_ptr, x_ptr, common_ptr->getMatStrain(),
1218 common_ptr->getMatCauchyStress(), vec_map,
1219 sym_map);
1220 };
1221
1222 auto get_tag_id_on_pmesh = [&](bool post_proc_skin) {
1223 int def_val_int = 0;
1224 Tag tag_mat;
1225 CHKERR mField.get_moab().tag_get_handle(
1226 "MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
1227 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
1228 auto meshset_vec_ptr =
1229 mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
1230 std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()));
1231
1232 Range skin_ents;
1233 std::unique_ptr<Skinner> skin_ptr;
1234 if (post_proc_skin) {
1235 skin_ptr = std::make_unique<Skinner>(&mField.get_moab());
1236 auto boundary_meshset = simple->getBoundaryMeshSet();
1237 CHKERR mField.get_moab().get_entities_by_handle(boundary_meshset,
1238 skin_ents, true);
1239 }
1240
1241 for (auto m : meshset_vec_ptr) {
1242 Range ents_3d;
1243 CHKERR mField.get_moab().get_entities_by_handle(m->getMeshset(), ents_3d,
1244 true);
1245 int const id = m->getMeshsetId();
1246 ents_3d = ents_3d.subset_by_dimension(SPACE_DIM);
1247 if (post_proc_skin) {
1248 Range skin_faces;
1249 CHKERR skin_ptr->find_skin(0, ents_3d, false, skin_faces);
1250 ents_3d = intersect(skin_ents, skin_faces);
1251 }
1252 CHKERR mField.get_moab().tag_clear_data(tag_mat, ents_3d, &id);
1253 }
1254
1255 return tag_mat;
1256 };
1257
1258 auto post_proc_domain = [&](auto post_proc_mesh) {
1259 auto post_proc_fe =
1260 boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
1262
1263 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
1264 calculate_stress_ops(post_proc_fe->getOpPtrVector());
1265
1266 post_proc_fe->getOpPtrVector().push_back(
1267
1268 new OpPPMap(
1269
1270 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1271
1272 {},
1273
1274 vec_map,
1275
1276 {},
1277
1278 sym_map
1279
1280 )
1281
1282 );
1283
1284 post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(false)});
1285 return post_proc_fe;
1286 };
1287
1288 auto post_proc_boundary = [&](auto post_proc_mesh) {
1289 auto post_proc_fe =
1290 boost::make_shared<PostProcEleBdy>(mField, post_proc_mesh);
1292 post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
1293 auto op_loop_side =
1294 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
1295 // push ops to side element, through op_loop_side operator
1296 auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
1297 calculate_stress_ops(op_loop_side->getOpPtrVector());
1298 post_proc_fe->getOpPtrVector().push_back(op_loop_side);
1299 auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
1300 post_proc_fe->getOpPtrVector().push_back(
1301 new ElasticExample::OpCalculateTraction(mat_stress_ptr,
1302 mat_traction_ptr));
1303 vec_map["T"] = mat_traction_ptr;
1304
1306
1307 post_proc_fe->getOpPtrVector().push_back(
1308
1309 new OpPPMap(
1310
1311 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
1312
1313 {},
1314
1315 vec_map,
1316
1317 {},
1318
1319 sym_map
1320
1321 )
1322
1323 );
1324
1325 post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(true)});
1326 return post_proc_fe;
1327 };
1328
1329 PetscBool post_proc_skin_only = PETSC_FALSE;
1330 if (SPACE_DIM == 3) {
1331 post_proc_skin_only = PETSC_TRUE;
1332 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin_only",
1333 &post_proc_skin_only, PETSC_NULLPTR);
1334 }
1335 if (post_proc_skin_only == PETSC_FALSE) {
1336 pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
1337 } else {
1338 pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
1339 }
1340
1342 post_proc_begin->getFEMethod());
1343 CHKERR pip->loopFiniteElements();
1345 post_proc_end->getFEMethod());
1346
1347 CHKERR post_proc_end->writeFile("out_elastic_" + std::to_string(iter) +
1348 ".h5m");
1350}
PetscErrorCode DMoFEMPreProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:536
virtual moab::Interface & get_moab()=0
Interface for managing meshsets containing materials and boundary conditions.
Element used to execute operators on side of the element.

◆ pushOperators()

MoFEMErrorCode Example::pushOperators ( )
private

[Set density distribution]

[Push operators to pipeline]

Definition at line 188 of file integration.cpp.

188 {
191
192 // Push an operator which calculates values of density at integration points
193 pipeline_mng->getOpDomainRhsPipeline().push_back(
195 "rho", commonDataPtr->getRhoAtIntegrationPtsPtr()));
196
197 // Push an operator to pipeline to calculate zero moment of inertia (mass)
198 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpZero(commonDataPtr));
199
200 // Push an operator to the pipeline to calculate first moment of inertaia
201 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpFirst(commonDataPtr));
202
203 // Push an operator to the pipeline to calculate second moment of inertaia
204 pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSecond(commonDataPtr));
205
206 // Set integration rule. Integration rule is equal to the polynomial order of
207 // the density field plus 2, since under the integral of the second moment of
208 // inertia term x*x is present
209 auto integration_rule = [](int, int, int p_data) { return p_data + 2; };
210
211 // Add integration rule to the element
214}
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)

◆ readMesh() [1/11]

MoFEMErrorCode Example::readMesh ( )
private

[Run problem]

[Run programme]

[run problem]

[Read mesh]

Examples
dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, mofem/tutorials/vec-7/adjoint.cpp, phase.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 462 of file dynamic_first_order_con_law.cpp.

462 {
465
467 CHKERR simple->loadFile();
469}
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180

◆ readMesh() [2/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [3/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [4/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [5/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [6/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [7/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [8/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [9/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [10/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [11/11]

MoFEMErrorCode Example::readMesh ( )
private

◆ rhsSource()

double Example::rhsSource ( const double x,
const double y,
const double  )
staticprivate
Examples
phase.cpp.

Definition at line 150 of file phase.cpp.

150 {
151 const auto idx = getCoordsInImage(x, y);
152
153 double v = 0;
154 for (auto w = 0; w != window_savitzky_golay; ++w) {
155 const auto i = focalIndex - (window_savitzky_golay - 1) / 2 + w;
156 const auto &intensity = iI[i];
157 v += intensity(idx.first, idx.second) * savitzkyGolayWeights[w];
158 }
159 v = static_cast<double>(v) / savitzkyGolayNormalisation;
160
161 const auto dz = rZ[focalIndex + 1] - rZ[focalIndex - 1];
162 return -k * v / dz;
163}
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'k', 3 > k
static int window_savitzky_golay
Definition phase.cpp:77
static int savitzkyGolayNormalisation
Definition phase.cpp:113
static std::vector< double > rZ
Definition phase.cpp:108
static const int * savitzkyGolayWeights
Definition phase.cpp:114

◆ runProblem() [1/14]

MoFEMErrorCode Example::runProblem ( )

[Run problem]

[Create common data]

[Run programme]

[Operator]

[run problem]

[Run all]

[Run problem]

Examples
dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, mofem/tutorials/vec-7/adjoint.cpp, phase.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 256 of file plastic.cpp.

256 {
260 CHKERR bC();
261 CHKERR OPs();
262 PetscBool test_ops = PETSC_FALSE;
263 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_operators", &test_ops,
264 PETSC_NULLPTR);
265 if (test_ops == PETSC_FALSE) {
266 CHKERR tsSolve();
267 } else {
269 }
271}
MoFEMErrorCode testOperators()
[Solve]
Definition plastic.cpp:1477
MoFEMErrorCode tsSolve()
Definition plastic.cpp:834
MoFEMErrorCode createCommonData()
[Set up problem]
Definition plastic.cpp:480
MoFEMErrorCode OPs()
[Boundary condition]
Definition plastic.cpp:650
MoFEMErrorCode setupProblem()
[Run problem]
Definition plastic.cpp:275
MoFEMErrorCode bC()
[Create common data]
Definition plastic.cpp:606

◆ runProblem() [2/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [3/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [4/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [5/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [6/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [7/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [8/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [9/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [10/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [11/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [12/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [13/14]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [14/14]

MoFEMErrorCode Example::runProblem ( )

◆ setFieldValues()

MoFEMErrorCode Example::setFieldValues ( )
private

[Create common data]

[Set density distribution]

Definition at line 172 of file integration.cpp.

172 {
174 auto set_density = [&](VectorAdaptor &&field_data, double *xcoord,
175 double *ycoord, double *zcoord) {
177 field_data[0] = 1;
179 };
180 FieldBlas *field_blas;
181 CHKERR mField.getInterface(field_blas);
182 CHKERR field_blas->setVertexDofs(set_density, "rho");
184}
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition Types.hpp:115
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.

◆ setIntegrationRules() [1/3]

MoFEMErrorCode Example::setIntegrationRules ( )
private

[Set up problem]

[Set integration rule]

Examples
heat_method.cpp, and plot_base.cpp.

Definition at line 229 of file plot_base.cpp.

229 {
232}

◆ setIntegrationRules() [2/3]

MoFEMErrorCode Example::setIntegrationRules ( )
private

◆ setIntegrationRules() [3/3]

MoFEMErrorCode Example::setIntegrationRules ( )
private

◆ setUp()

MoFEMErrorCode Example::setUp ( )
private

[Run all]

[Set up problem]

Definition at line 137 of file integration.cpp.

137 {
141 CHKERR simple->loadFile();
142 // Add field
143 CHKERR simple->addDomainField("rho", H1, AINSWORTH_LEGENDRE_BASE, 1);
144 constexpr int order = 1;
145 CHKERR simple->setFieldOrder("rho", order);
146 CHKERR simple->setUp();
148}
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base nme:nme847.
Definition definitions.h:60

◆ setupAdJoint()

MoFEMErrorCode Example::setupAdJoint ( )
private

[Set up problem]

[Setup adjoint]

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 439 of file adjoint.cpp.

439 {
442
443 // ad adjoint dm
444
445 auto create_adjoint_dm = [&]() {
446 auto adjoint_dm = createDM(mField.get_comm(), "DMMOFEM");
447
448 auto add_field = [&]() {
450 CHKERR mField.add_field("ADJOINT_FIELD", H1, base, SPACE_DIM);
452 "ADJOINT_FIELD");
453 for (auto tt = MBEDGE; tt <= moab::CN::TypeDimensionMap[SPACE_DIM].second;
454 ++tt)
455 CHKERR mField.set_field_order(simple->getMeshset(), tt, "ADJOINT_FIELD",
456 fieldOrder);
457 CHKERR mField.set_field_order(simple->getMeshset(), MBVERTEX,
458 "ADJOINT_FIELD", 1);
461 };
462
463 auto add_adjoint_fe_impl = [&]() {
465 CHKERR mField.add_finite_element("ADJOINT_DOMAIN_FE");
467 "ADJOINT_FIELD");
469 "ADJOINT_FIELD");
471 "ADJOINT_FIELD");
473 "GEOMETRY");
475 simple->getMeshset(), SPACE_DIM, "ADJOINT_DOMAIN_FE");
476 CHKERR mField.build_finite_elements("ADJOINT_DOMAIN_FE");
477
478 CHKERR mField.add_finite_element("ADJOINT_BOUNDARY_FE");
480 "ADJOINT_FIELD");
482 "ADJOINT_FIELD");
484 "ADJOINT_FIELD");
486 "GEOMETRY");
487
488 auto block_name = "OPTIMISE";
489 auto mesh_mng = mField.getInterface<MeshsetsManager>();
490 auto bcs = mesh_mng->getCubitMeshsetPtr(
491
492 std::regex((boost::format("%s(.*)") % block_name).str())
493
494 );
495
496 for (auto bc : bcs) {
498 bc->getMeshset(), SPACE_DIM - 1, "ADJOINT_BOUNDARY_FE");
499 }
500
501 CHKERR mField.build_finite_elements("ADJOINT_BOUNDARY_FE");
502
503 CHKERR mField.build_adjacencies(simple->getBitRefLevel(),
504 simple->getBitRefLevelMask());
505
507 };
508
509 auto set_adjoint_dm_imp = [&]() {
511 CHKERR DMMoFEMCreateMoFEM(adjoint_dm, &mField, "ADJOINT",
512 simple->getBitRefLevel(),
513 simple->getBitRefLevelMask());
514 CHKERR DMMoFEMSetDestroyProblem(adjoint_dm, PETSC_TRUE);
515 CHKERR DMSetFromOptions(adjoint_dm);
516 CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_DOMAIN_FE");
517 CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_BOUNDARY_FE");
518 CHKERR DMMoFEMSetSquareProblem(adjoint_dm, PETSC_TRUE);
519 CHKERR DMMoFEMSetIsPartitioned(adjoint_dm, PETSC_TRUE);
520 mField.getInterface<ProblemsManager>()->buildProblemFromFields =
521 PETSC_TRUE;
522 CHKERR DMSetUp(adjoint_dm);
523 mField.getInterface<ProblemsManager>()->buildProblemFromFields =
524 PETSC_FALSE;
526 };
527
528 CHK_THROW_MESSAGE(add_field(), "add adjoint field");
529 CHK_THROW_MESSAGE(add_adjoint_fe_impl(), "add adjoint fe");
530 CHK_THROW_MESSAGE(set_adjoint_dm_imp(), "set adjoint dm");
531
532 return adjoint_dm;
533 };
534
535 adjointDM = create_adjoint_dm();
536
538}
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
PetscErrorCode DMMoFEMSetIsPartitioned(DM dm, PetscBool is_partitioned)
Definition DMMoFEM.cpp:1113
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition DMMoFEM.cpp:488
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition DMMoFEM.cpp:450
PetscErrorCode DMMoFEMCreateMoFEM(DM dm, MoFEM::Interface *m_field_ptr, const char problem_name[], const MoFEM::BitRefLevel bit_level, const MoFEM::BitRefLevel bit_mask=MoFEM::BitRefLevel().set())
Must be called by user to set MoFEM data structures.
Definition DMMoFEM.cpp:114
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string name, const bool recursive=true)=0
add entities to finite element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
virtual MoFEMErrorCode add_ents_to_field_by_dim(const Range &ents, const int dim, const std::string &name, int verb=DEFAULT_VERBOSITY)=0
Add entities to field meshset.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
PetscErrorCode DMMoFEMSetDestroyProblem(DM dm, PetscBool destroy_problem)
Definition DMMoFEM.cpp:434
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
FieldApproximationBase base
Definition plot_base.cpp:68
int fieldOrder
Definition adjoint.cpp:153
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
virtual MoFEMErrorCode add_field(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_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
Problem manager is used to build and partition problems.

◆ setupProblem() [1/13]

MoFEMErrorCode Example::setupProblem ( )
private

[Run problem]

[Read mesh]

[Set up problem]

Examples
dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, mofem/tutorials/vec-7/adjoint.cpp, phase.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 275 of file plastic.cpp.

275 {
278
279 Range domain_ents;
280 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
281 true);
282 auto get_ents_by_dim = [&](const auto dim) {
283 if (dim == SPACE_DIM) {
284 return domain_ents;
285 } else {
286 Range ents;
287 if (dim == 0)
288 CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
289 else
290 CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
291 return ents;
292 }
293 };
294
295 auto get_base = [&]() {
296 auto domain_ents = get_ents_by_dim(SPACE_DIM);
297 if (domain_ents.empty())
298 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
299 const auto type = type_from_handle(domain_ents[0]);
300 switch (type) {
301 case MBQUAD:
303 case MBHEX:
305 case MBTRI:
307 case MBTET:
309 default:
310 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
311 }
312 return NOBASE;
313 };
314
315 const auto base = get_base();
316 MOFEM_LOG("PLASTICITY", Sev::inform)
317 << "Base " << ApproximationBaseNames[base];
318
319 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
320 CHKERR simple->addDomainField("EP", L2, base, size_symm);
321 CHKERR simple->addDomainField("TAU", L2, base, 1);
322 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
323
324 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
325
326 PetscBool order_edge = PETSC_FALSE;
327 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-order_edge", &order_edge,
328 PETSC_NULLPTR);
329 PetscBool order_face = PETSC_FALSE;
330 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-order_face", &order_face,
331 PETSC_NULLPTR);
332 PetscBool order_volume = PETSC_FALSE;
333 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-order_volume", &order_volume,
334 PETSC_NULLPTR);
335
336 if (order_edge || order_face || order_volume) {
337
338 MOFEM_LOG("PLASTICITY", Sev::inform) << "Order edge " << order_edge
339 ? "true"
340 : "false";
341 MOFEM_LOG("PLASTICITY", Sev::inform) << "Order face " << order_face
342 ? "true"
343 : "false";
344 MOFEM_LOG("PLASTICITY", Sev::inform) << "Order volume " << order_volume
345 ? "true"
346 : "false";
347
348 auto ents = get_ents_by_dim(0);
349 if (order_edge)
350 ents.merge(get_ents_by_dim(1));
351 if (order_face)
352 ents.merge(get_ents_by_dim(2));
353 if (order_volume)
354 ents.merge(get_ents_by_dim(3));
355 CHKERR simple->setFieldOrder("U", order, &ents);
356 } else {
357 CHKERR simple->setFieldOrder("U", order);
358 }
359 CHKERR simple->setFieldOrder("EP", ep_order);
360 CHKERR simple->setFieldOrder("TAU", tau_order);
361
362 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
363
364#ifdef ADD_CONTACT
365 CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
366 SPACE_DIM);
367 CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
368 SPACE_DIM);
369
370 auto get_skin = [&]() {
371 Range body_ents;
372 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
373 Skinner skin(&mField.get_moab());
374 Range skin_ents;
375 CHKERR skin.find_skin(0, body_ents, false, skin_ents);
376 return skin_ents;
377 };
378
379 auto filter_blocks = [&](auto skin) {
380 bool is_contact_block = true;
381 Range contact_range;
382 for (auto m :
384
385 (boost::format("%s(.*)") % "CONTACT").str()
386
387 ))
388
389 ) {
390 is_contact_block =
391 true; ///< blocs interation is collective, so that is set irrespective
392 ///< if there are entities in given rank or not in the block
393 MOFEM_LOG("CONTACT", Sev::inform)
394 << "Find contact block set: " << m->getName();
395 auto meshset = m->getMeshset();
396 Range contact_meshset_range;
397 CHKERR mField.get_moab().get_entities_by_dimension(
398 meshset, SPACE_DIM - 1, contact_meshset_range, true);
399
400 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
401 contact_meshset_range);
402 contact_range.merge(contact_meshset_range);
403 }
404 if (is_contact_block) {
405 MOFEM_LOG("SYNC", Sev::inform)
406 << "Nb entities in contact surface: " << contact_range.size();
408 skin = intersect(skin, contact_range);
409 }
410 return skin;
411 };
412
413 auto filter_true_skin = [&](auto skin) {
414 Range boundary_ents;
415 ParallelComm *pcomm =
416 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
417 CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
418 PSTATUS_NOT, -1, &boundary_ents);
419 return boundary_ents;
420 };
421
422 auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
423 CHKERR simple->setFieldOrder("SIGMA", 0);
424 CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
425#endif
426
427 CHKERR simple->setUp();
428 CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
429
430 auto project_ho_geometry = [&]() {
431 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
432 return mField.loop_dofs("GEOMETRY", ent_method);
433 };
434 PetscBool project_geometry = PETSC_TRUE;
435 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-project_geometry",
436 &project_geometry, PETSC_NULLPTR);
437 if (project_geometry) {
438 CHKERR project_ho_geometry();
439 }
440
441 auto get_volume = [&]() {
442 using VolOp = DomainEle::UserDataOperator;
443 auto *op_ptr = new VolOp(NOSPACE, VolOp::OPSPACE);
444 std::array<double, 2> volume_and_count;
445 op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
446 EntityType type,
449 auto op_ptr = static_cast<VolOp *>(base_op_ptr);
450 volume_and_count[VOL] += op_ptr->getMeasure();
451 volume_and_count[COUNT] += 1;
452 // in necessary at integration over Gauss points.
454 };
455 volume_and_count = {0, 0};
456 auto fe = boost::make_shared<DomainEle>(mField);
457 fe->getOpPtrVector().push_back(op_ptr);
458
459 auto dm = simple->getDM();
461 DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), fe),
462 "cac volume");
463 std::array<double, 2> tot_volume_and_count;
464 MPI_Allreduce(volume_and_count.data(), tot_volume_and_count.data(),
465 volume_and_count.size(), MPI_DOUBLE, MPI_SUM,
466 mField.get_comm());
467 return tot_volume_and_count;
468 };
469
470 meshVolumeAndCount = get_volume();
471 MOFEM_LOG("PLASTICITY", Sev::inform)
472 << "Mesh volume " << meshVolumeAndCount[VOL] << " nb. of elements "
474
476}
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
@ NOBASE
Definition definitions.h:59
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
@ L2
field with C-1 continuity
Definition definitions.h:88
@ NOSPACE
Definition definitions.h:83
#define MYPCOMM_INDEX
default communicator number PCOMM
@ MOFEM_NOT_FOUND
Definition definitions.h:33
static const char *const ApproximationBaseNames[]
Definition definitions.h:72
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
auto type_from_handle(const EntityHandle h)
get type from entity handle
constexpr auto size_symm
Definition plastic.cpp:42
constexpr FieldSpace CONTACT_SPACE
Definition plastic.cpp:52
static std::array< double, 2 > meshVolumeAndCount
Definition plastic.cpp:223
Managing BitRefLevels.
base operator to do operations at Gauss Pt. level
Data on single entity (This is passed as argument to DataOperator::doWork)
double getMeasure() const
get measure of element
@ OPSPACE
operator do Work is execute on space data
Projection of edge entities with one mid-node on hierarchical basis.
VolEle::UserDataOperator VolOp

◆ setupProblem() [2/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [3/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [4/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [5/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [6/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [7/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [8/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [9/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [10/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [11/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [12/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [13/13]

MoFEMErrorCode Example::setupProblem ( )
private

◆ solveElastic()

MoFEMErrorCode Example::solveElastic ( )
private

[Push operators to pipeline]

[Solve]

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 1036 of file adjoint.cpp.

1036 {
1038 auto simple = mField.getInterface<Simple>();
1039 auto dm = simple->getDM();
1040 auto f = createDMVector(dm);
1041 auto d = vectorDuplicate(f);
1042 CHKERR VecZeroEntries(d);
1043 CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
1044
1045 auto set_essential_bc = [&]() {
1047 // This is low level pushing finite elements (pipelines) to solver
1048
1049 auto ksp_ctx_ptr = getDMKspCtx(dm);
1050 auto pre_proc_rhs = boost::make_shared<FEMethod>();
1051 auto post_proc_rhs = boost::make_shared<FEMethod>();
1052 auto post_proc_lhs = boost::make_shared<FEMethod>();
1053
1054 auto get_pre_proc_hook = [&]() {
1056 {});
1057 };
1058 pre_proc_rhs->preProcessHook = get_pre_proc_hook();
1059
1060 auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
1062
1064 post_proc_rhs, 1.)();
1066 };
1067
1068 auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
1070
1072 post_proc_lhs, 1.)();
1074 };
1075
1076 post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
1077 post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
1078
1079 ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
1080 ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
1081 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
1083 };
1084
1085 auto setup_and_solve = [&](auto solver) {
1087 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1088 MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp";
1089 CHKERR KSPSetUp(solver);
1090 MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp <= Done";
1091 MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve";
1092 CHKERR KSPSolve(solver, f, d);
1093 MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve <= Done";
1095 };
1096
1097 MOFEM_LOG_CHANNEL("TIMER");
1098 MOFEM_LOG_TAG("TIMER", "timer");
1099
1100 CHKERR set_essential_bc();
1101 CHKERR setup_and_solve(kspElastic);
1102
1103 CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
1104 CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
1105 CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
1106
1107 auto evaluate_field_at_the_point = [&]() {
1109
1110 int coords_dim = 3;
1111 std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
1112 PetscBool do_eval_field = PETSC_FALSE;
1113 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
1114 field_eval_coords.data(), &coords_dim,
1115 &do_eval_field);
1116
1117 if (do_eval_field) {
1118
1119 vectorFieldPtr = boost::make_shared<MatrixDouble>();
1120 auto field_eval_data =
1121 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
1122
1124 ->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
1125
1126 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1127 auto no_rule = [](int, int, int) { return -1; };
1128 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1129 field_eval_fe_ptr->getRuleHook = no_rule;
1130
1131 field_eval_fe_ptr->getOpPtrVector().push_back(
1133
1135 ->evalFEAtThePoint<SPACE_DIM>(
1136 field_eval_coords.data(), 1e-12, simple->getProblemName(),
1137 simple->getDomainFEName(), field_eval_data,
1139 QUIET);
1140
1141 if (vectorFieldPtr->size1()) {
1143 if constexpr (SPACE_DIM == 2)
1144 MOFEM_LOG("FieldEvaluator", Sev::inform)
1145 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
1146 else
1147 MOFEM_LOG("FieldEvaluator", Sev::inform)
1148 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
1149 << " U_Z: " << t_disp(2);
1150 }
1151
1153 }
1155 };
1156
1157 CHKERR evaluate_field_at_the_point();
1158
1160}
@ QUIET
@ MF_EXIST
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
PetscErrorCode PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
auto getDMKspCtx(DM dm)
Get KSP context data structure used by DM.
Definition DMMoFEM.hpp:1116
PetscBool do_eval_field
Evaluate field.
Definition plastic.cpp:119
boost::shared_ptr< MatrixDouble > vectorFieldPtr
Definition elastic.cpp:87
Class (Function) to enforce essential constrains on the left hand side diagonal.
Definition Essential.hpp:33
Field evaluator interface.
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:800

◆ solveSystem() [1/10]

MoFEMErrorCode Example::solveSystem ( )
private

[Solve]

[Push operators to pipeline]

[Solve]

< Mass matrix

< Linear solver

Examples
dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, phase.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 877 of file dynamic_first_order_con_law.cpp.

877 {
879 auto *simple = mField.getInterface<Simple>();
880 auto *pipeline_mng = mField.getInterface<PipelineManager>();
881
882 auto dm = simple->getDM();
883
884 auto calculate_stress_ops = [&](auto &pip) {
886
887 auto v_ptr = boost::make_shared<MatrixDouble>();
888 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("V", v_ptr));
889 auto X_ptr = boost::make_shared<MatrixDouble>();
890 pip.push_back(
891 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
892
893 auto x_ptr = boost::make_shared<MatrixDouble>();
894 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("x_1", x_ptr));
895
896 // Calculate unknown F
897 auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
899 "F", mat_H_tensor_ptr));
900
901 auto u_ptr = boost::make_shared<MatrixDouble>();
902 pip.push_back(new OpCalculateDisplacement<SPACE_DIM>(x_ptr, X_ptr, u_ptr));
903 // Calculate P
904
905 auto mat_F_ptr = boost::make_shared<MatrixDouble>();
907 mat_F_ptr, mat_H_tensor_ptr));
908
909 PetscBool is_linear_elasticity = PETSC_TRUE;
910 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_linear_elasticity",
911 &is_linear_elasticity, PETSC_NULLPTR);
912
913 auto mat_P_ptr = boost::make_shared<MatrixDouble>();
914 if (is_linear_elasticity) {
917 mat_F_ptr));
918 } else {
919 auto inv_F = boost::make_shared<MatrixDouble>();
920 auto det_ptr = boost::make_shared<VectorDouble>();
921
922 pip.push_back(new OpInvertMatrix<SPACE_DIM>(mat_F_ptr, det_ptr, inv_F));
923
926 mat_F_ptr, inv_F, det_ptr));
927 }
928
929 auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
931 "V", mat_v_grad_ptr));
932
933 return boost::make_tuple(v_ptr, X_ptr, x_ptr, mat_P_ptr, mat_F_ptr, u_ptr);
934 };
935
936 auto post_proc_boundary = [&]() {
937 auto boundary_post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
938
940 boundary_post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
941 auto op_loop_side =
942 new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
943 // push ops to side element, through op_loop_side operator
944 auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
945 boundary_mat_F_ptr, boundary_u_ptr] =
946 calculate_stress_ops(op_loop_side->getOpPtrVector());
947 boundary_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
948
950
951 boundary_post_proc_fe->getOpPtrVector().push_back(
952
953 new OpPPMap(
954
955 boundary_post_proc_fe->getPostProcMesh(),
956 boundary_post_proc_fe->getMapGaussPts(),
957
959
960 OpPPMap::DataMapMat{{"V", boundary_v_ptr},
961 {"GEOMETRY", boundary_X_ptr},
962 {"x", boundary_x_ptr},
963 {"U", boundary_u_ptr}},
964
965 OpPPMap::DataMapMat{{"FIRST_PIOLA", boundary_mat_P_ptr},
966 {"F", boundary_mat_F_ptr}},
967
969
970 )
971
972 );
973 return boundary_post_proc_fe;
974 };
975
976 // Add monitor to time solver
977
978 double rho = 1.;
979 CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-density", &rho, PETSC_NULLPTR);
980 auto get_rho = [rho](const double, const double, const double) {
981 return rho;
982 };
983
984 SmartPetscObj<Mat> M; ///< Mass matrix
985 SmartPetscObj<KSP> ksp; ///< Linear solver
986
987 auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
988 tsPrePostProc = ts_pre_post_proc;
989
991 CHKERR MatZeroEntries(M);
992
993 boost::shared_ptr<DomainEle> vol_mass_ele(new DomainEle(mField));
994
995 vol_mass_ele->B = M;
996
997 auto integration_rule = [](int, int, int approx_order) {
998 return 2 * approx_order;
999 };
1000
1001 vol_mass_ele->getRuleHook = integration_rule;
1002
1003 vol_mass_ele->getOpPtrVector().push_back(new OpMassV("V", "V", get_rho));
1004 vol_mass_ele->getOpPtrVector().push_back(new OpMassF("F", "F"));
1005
1006 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), vol_mass_ele);
1007 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
1008 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
1009
1010 auto lumpVec = createDMVector(simple->getDM());
1011 CHKERR MatGetRowSum(M, lumpVec);
1012
1013 CHKERR MatZeroEntries(M);
1014 CHKERR MatDiagonalSet(M, lumpVec, INSERT_VALUES);
1015
1016 // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
1017 // M^-1G(t,u)
1018 ksp = createKSP(mField.get_comm());
1019 CHKERR KSPSetOperators(ksp, M, M);
1020 CHKERR KSPSetFromOptions(ksp);
1021 CHKERR KSPSetUp(ksp);
1022
1023 auto solve_boundary_for_g = [&]() {
1025 if (*(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch)) {
1026
1027 CHKERR VecGhostUpdateBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
1028 ADD_VALUES, SCATTER_REVERSE);
1029 CHKERR VecGhostUpdateEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
1030 ADD_VALUES, SCATTER_REVERSE);
1031 CHKERR VecAssemblyBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1032 CHKERR VecAssemblyEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1033 *(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch) = false;
1034
1035 auto D =
1036 smartVectorDuplicate(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1037 CHKERR KSPSolve(ksp, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F, D);
1038 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1039 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1040 CHKERR VecCopy(D, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
1041 }
1042
1044 };
1045
1046 pipeline_mng->getBoundaryExplicitRhsFE()->postProcessHook =
1047 solve_boundary_for_g;
1048
1050 ts = pipeline_mng->createTSEX(dm);
1051
1052 // Field eval
1053 PetscBool field_eval_flag = PETSC_TRUE;
1054 boost::shared_ptr<MatrixDouble> velocity_field_ptr;
1055 boost::shared_ptr<MatrixDouble> geometry_field_ptr;
1056 boost::shared_ptr<MatrixDouble> spatial_position_field_ptr;
1057 boost::shared_ptr<SetPtsData> field_eval_data;
1058
1059 std::array<double, 3> field_eval_coords = {0.5, 0.5, 5.};
1060 int dim = 3;
1061 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
1062 field_eval_coords.data(), &dim,
1063 &field_eval_flag);
1064
1065 if (field_eval_flag) {
1066 field_eval_data =
1067 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
1068 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
1069 field_eval_data, simple->getDomainFEName());
1070
1071 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1072
1073 auto no_rule = [](int, int, int) { return -1; };
1074
1075 auto fe_ptr = field_eval_data->feMethodPtr.lock();
1076 fe_ptr->getRuleHook = no_rule;
1077 velocity_field_ptr = boost::make_shared<MatrixDouble>();
1078 geometry_field_ptr = boost::make_shared<MatrixDouble>();
1079 spatial_position_field_ptr = boost::make_shared<MatrixDouble>();
1080 fe_ptr->getOpPtrVector().push_back(
1081 new OpCalculateVectorFieldValues<SPACE_DIM>("V", velocity_field_ptr));
1082 fe_ptr->getOpPtrVector().push_back(
1084 geometry_field_ptr));
1085 fe_ptr->getOpPtrVector().push_back(
1087 "x_2", spatial_position_field_ptr));
1088 }
1089
1090 auto post_proc_domain = [&]() {
1091 auto post_proc_fe_vol = boost::make_shared<PostProcEle>(mField);
1092
1094
1095 auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
1096 boundary_mat_F_ptr, boundary_u_ptr] =
1097 calculate_stress_ops(post_proc_fe_vol->getOpPtrVector());
1098
1099 post_proc_fe_vol->getOpPtrVector().push_back(
1100
1101 new OpPPMap(
1102
1103 post_proc_fe_vol->getPostProcMesh(),
1104 post_proc_fe_vol->getMapGaussPts(),
1105
1106 {},
1107
1108 {{"V", boundary_v_ptr},
1109 {"GEOMETRY", boundary_X_ptr},
1110 {"x", boundary_x_ptr},
1111 {"U", boundary_u_ptr}},
1112
1113 {{"FIRST_PIOLA", boundary_mat_P_ptr}, {"F", boundary_mat_F_ptr}},
1114
1115 {}
1116
1117 )
1118
1119 );
1120 return post_proc_fe_vol;
1121 };
1122
1123 boost::shared_ptr<FEMethod> null_fe;
1124 auto monitor_ptr = boost::make_shared<Monitor>(
1125 SmartPetscObj<DM>(dm, true), mField, post_proc_domain(),
1126 post_proc_boundary(), velocity_field_ptr, spatial_position_field_ptr,
1127 geometry_field_ptr, field_eval_coords, field_eval_data);
1128
1129 CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
1130 null_fe, monitor_ptr);
1131
1132 double ftime = 1;
1133 // CHKERR TSSetMaxTime(ts, ftime);
1134 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
1135
1136 auto T = createDMVector(simple->getDM());
1137 CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
1138 SCATTER_FORWARD);
1139 CHKERR TSSetSolution(ts, T);
1140 CHKERR TSSetFromOptions(ts);
1141
1142 CHKERR TSSetPostStage(ts, TSPrePostProc::tsPostStage);
1143 CHKERR TSSetPostStep(ts, TSPrePostProc::tsPostStep);
1144 CHKERR TSSetPreStep(ts, TSPrePostProc::tsPreStep);
1145
1146 boost::shared_ptr<ForcesAndSourcesCore> null;
1147
1148 if (auto ptr = tsPrePostProc.lock()) {
1149 ptr->fsRawPtr = this;
1150 CHKERR TSSetUp(ts);
1151 CHKERR TSSolve(ts, NULL);
1152 CHKERR TSGetTime(ts, &ftime);
1153 }
1154
1156}
PipelineManager::ElementsAndOpsByDim< SPACE_DIM >::DomainEle DomainEle
Definition adjoint.cpp:32
static boost::weak_ptr< TSPrePostProc > tsPrePostProc
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMassV
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM *SPACE_DIM > OpMassF
PetscErrorCode DMCreateMatrix_MoFEM(DM dm, Mat *M)
Definition DMMoFEM.cpp:1187
double D
auto createKSP(MPI_Comm comm)
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.
Definition DMMoFEM.cpp:1046
DEPRECATED SmartPetscObj< Vec > smartVectorDuplicate(Vec vec)
SmartPetscObj< Mat > M
intrusive_ptr for managing petsc objects
static MoFEMErrorCode tsPostStep(TS ts)
static MoFEMErrorCode tsPreStep(TS ts)
static MoFEMErrorCode tsPostStage(TS ts, PetscReal stagetime, PetscInt stageindex, Vec *Y)
[Boundary condition]

◆ solveSystem() [2/10]

MoFEMErrorCode Example::solveSystem ( )
private

◆ solveSystem() [3/10]

MoFEMErrorCode Example::solveSystem ( )
private

◆ solveSystem() [4/10]

MoFEMErrorCode Example::solveSystem ( )
private

◆ solveSystem() [5/10]

MoFEMErrorCode Example::solveSystem ( )
private

◆ solveSystem() [6/10]

MoFEMErrorCode Example::solveSystem ( )
private

◆ solveSystem() [7/10]

MoFEMErrorCode Example::solveSystem ( )
private

◆ solveSystem() [8/10]

MoFEMErrorCode Example::solveSystem ( )
private

◆ solveSystem() [9/10]

MoFEMErrorCode Example::solveSystem ( )
private

◆ solveSystem() [10/10]

MoFEMErrorCode Example::solveSystem ( )
private

◆ testOperators()

MoFEMErrorCode Example::testOperators ( )
private

[Solve]

[TestOperators]

Examples
plastic.cpp.

Definition at line 1477 of file plastic.cpp.

1477 {
1479
1480 // get operators tester
1481 auto simple = mField.getInterface<Simple>();
1482 auto opt = mField.getInterface<OperatorsTester>(); // get interface to
1483 // OperatorsTester
1484 auto pip = mField.getInterface<PipelineManager>(); // get interface to
1485 // pipeline manager
1486
1487 constexpr double eps = 1e-9;
1488
1489 auto x = opt->setRandomFields(simple->getDM(), {
1490
1491 {"U", {-1e-6, 1e-6}},
1492
1493 {"EP", {-1e-6, 1e-6}},
1494
1495 {"TAU", {0, 1e-4}}
1496
1497 });
1498
1499 auto dot_x_plastic_active =
1500 opt->setRandomFields(simple->getDM(), {
1501
1502 {"U", {-1, 1}},
1503
1504 {"EP", {-1, 1}},
1505
1506 {"TAU", {0.1, 1}}
1507
1508 });
1509 auto diff_x_plastic_active =
1510 opt->setRandomFields(simple->getDM(), {
1511
1512 {"U", {-1, 1}},
1513
1514 {"EP", {-1, 1}},
1515
1516 {"TAU", {-1, 1}}
1517
1518 });
1519
1520 auto dot_x_elastic =
1521 opt->setRandomFields(simple->getDM(), {
1522
1523 {"U", {-1, 1}},
1524
1525 {"EP", {-1, 1}},
1526
1527 {"TAU", {-1, -0.1}}
1528
1529 });
1530 auto diff_x_elastic =
1531 opt->setRandomFields(simple->getDM(), {
1532
1533 {"U", {-1, 1}},
1534
1535 {"EP", {-1, 1}},
1536
1537 {"TAU", {-1, 1}}
1538
1539 });
1540
1541 auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline, auto rhs_pipeline,
1542 auto dot_x, auto diff_x) {
1544
1545 auto diff_res = opt->checkCentralFiniteDifference(
1546 simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
1547 SmartPetscObj<Vec>(), diff_x, 0, 0.5, eps);
1548
1549 // Calculate norm of difference between directional derivative calculated
1550 // from finite difference, and tangent matrix.
1551 double fnorm;
1552 CHKERR VecNorm(diff_res, NORM_2, &fnorm);
1553 MOFEM_LOG_C("PLASTICITY", Sev::inform,
1554 "Test consistency of tangent matrix %3.4e", fnorm);
1555
1556 constexpr double err = 1e-5;
1557 if (fnorm > err)
1558 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1559 "Norm of directional derivative too large err = %3.4e", fnorm);
1560
1562 };
1563
1564 MOFEM_LOG("PLASTICITY", Sev::inform) << "Elastic active";
1565 CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1566 pip->getDomainRhsFE(), dot_x_elastic, diff_x_elastic);
1567
1568 MOFEM_LOG("PLASTICITY", Sev::inform) << "Plastic active";
1569 CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
1570 pip->getDomainRhsFE(), dot_x_plastic_active,
1571 diff_x_plastic_active);
1572
1574};
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...

◆ topologyModes()

MoFEMErrorCode Example::topologyModes ( )
private

[Boundary condition]

[Adjoint modes]

Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 563 of file adjoint.cpp.

563 {
565
566 auto opt_ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
567 auto subset_dm_bdy = createDM(mField.get_comm(), "DMMOFEM");
568 CHKERR DMMoFEMSetSquareProblem(subset_dm_bdy, PETSC_TRUE);
569 CHKERR DMMoFEMCreateSubDM(subset_dm_bdy, adjointDM, "SUBSET_BDY");
570 CHKERR DMMoFEMAddElement(subset_dm_bdy, "ADJOINT_BOUNDARY_FE");
571 CHKERR DMMoFEMAddSubFieldRow(subset_dm_bdy, "ADJOINT_FIELD",
572 boost::make_shared<Range>(opt_ents));
573 CHKERR DMMoFEMAddSubFieldCol(subset_dm_bdy, "ADJOINT_FIELD",
574 boost::make_shared<Range>(opt_ents));
575 CHKERR DMSetUp(subset_dm_bdy);
576
577 auto subset_dm_domain = createDM(mField.get_comm(), "DMMOFEM");
578 CHKERR DMMoFEMSetSquareProblem(subset_dm_domain, PETSC_TRUE);
579 CHKERR DMMoFEMCreateSubDM(subset_dm_domain, adjointDM, "SUBSET_DOMAIN");
580 CHKERR DMMoFEMAddElement(subset_dm_domain, "ADJOINT_DOMAIN_FE");
581 CHKERR DMMoFEMAddSubFieldRow(subset_dm_domain, "ADJOINT_FIELD");
582 CHKERR DMMoFEMAddSubFieldCol(subset_dm_domain, "ADJOINT_FIELD");
583 CHKERR DMSetUp(subset_dm_domain);
584
585 // remove dofs on boundary of the domain
586 auto remove_dofs = [&]() {
588
589 std::array<Range, 3> remove_dim_ents;
590 remove_dim_ents[0] =
591 get_range_from_block(mField, "OPT_REMOVE_X", SPACE_DIM - 1);
592 remove_dim_ents[1] =
593 get_range_from_block(mField, "OPT_REMOVE_Y", SPACE_DIM - 1);
594 remove_dim_ents[2] =
595 get_range_from_block(mField, "OPT_REMOVE_Z", SPACE_DIM - 1);
596
597 for (int d = 0; d != 3; ++d) {
598 MOFEM_LOG("WORLD", Sev::inform)
599 << "Removing topology modes on block OPT_REMOVE_" << (char)('X' + d)
600 << " with " << remove_dim_ents[d].size() << " entities";
601 }
602
603 Range body_ents;
604 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents,
605 true);
606 auto skin = moab::Skinner(&mField.get_moab());
607 Range boundary_ents;
608 CHKERR skin.find_skin(0, body_ents, false, boundary_ents);
609 for (int d = 0; d != 3; ++d) {
610 boundary_ents = subtract(boundary_ents, remove_dim_ents[d]);
611 }
612 ParallelComm *pcomm =
613 ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
614 CHKERR pcomm->filter_pstatus(boundary_ents,
615 PSTATUS_SHARED | PSTATUS_MULTISHARED,
616 PSTATUS_NOT, -1, nullptr);
617 for (auto d = SPACE_DIM - 2; d >= 0; --d) {
618 if (d >= 0) {
619 Range ents;
620 CHKERR mField.get_moab().get_adjacencies(boundary_ents, d, false, ents,
621 moab::Interface::UNION);
622 boundary_ents.merge(ents);
623 } else {
624 Range verts;
625 CHKERR mField.get_moab().get_connectivity(boundary_ents, verts);
626 boundary_ents.merge(verts);
627 }
628 CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
629 boundary_ents);
630 }
631 boundary_ents.merge(opt_ents);
632 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
633 "SUBSET_DOMAIN", "ADJOINT_FIELD", boundary_ents);
634 for (int d = 0; d != 3; ++d) {
635 CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
636 "SUBSET_DOMAIN", "ADJOINT_FIELD", remove_dim_ents[d], d, d);
637 }
638
639 // #ifndef NDEBUG
640 if (mField.get_comm_rank() == 0) {
641 CHKERR save_range(mField.get_moab(), "topoMode_boundary_ents.vtk",
642 boundary_ents);
643 }
644 // #endif
645
647 };
648
649 CHKERR remove_dofs();
650
651 auto get_lhs_fe = [&]() {
652 auto fe_lhs = boost::make_shared<BoundaryEle>(mField);
653 fe_lhs->getRuleHook = [](int, int, int p_data) {
654 return 2 * p_data + p_data - 1;
655 };
656 auto &pip = fe_lhs->getOpPtrVector();
658 "GEOMETRY");
661 pip.push_back(new OpMass("ADJOINT_FIELD", "ADJOINT_FIELD",
662 [](double, double, double) { return 1.; }));
663 return fe_lhs;
664 };
665
666 auto get_rhs_fe = [&]() {
667 auto fe_rhs = boost::make_shared<BoundaryEle>(mField);
668 fe_rhs->getRuleHook = [](int, int, int p_data) {
669 return 2 * p_data + p_data - 1;
670 };
671 auto &pip = fe_rhs->getOpPtrVector();
673 "GEOMETRY");
674
675 return fe_rhs;
676 };
677
678 auto block_name = "OPTIMISE";
679 auto mesh_mng = mField.getInterface<MeshsetsManager>();
680 auto bcs = mesh_mng->getCubitMeshsetPtr(
681
682 std::regex((boost::format("%s(.*)") % block_name).str())
683
684 );
685
686 for (auto &v : modeVecs) {
687 v = createDMVector(subset_dm_bdy);
688 }
689
691 struct OpMode : public OP {
692 OpMode(const std::string name,
693 boost::shared_ptr<ObjectiveFunctionData> python_ptr, int id,
694 std::vector<SmartPetscObj<Vec>> mode_vecs,
695 std::vector<std::array<double, 3>> mode_centroids,
696 std::vector<std::array<double, 6>> mode_bboxes, int block_counter,
697 int mode_counter, boost::shared_ptr<Range> range = nullptr)
698 : OP(name, name, OP::OPROW, range), pythonPtr(python_ptr), iD(id),
699 modeVecs(mode_vecs), modeCentroids(mode_centroids),
700 modeBboxes(mode_bboxes), blockCounter(block_counter),
701 modeCounter(mode_counter) {}
702
703 MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
705
706 if (OP::entsPtr) {
707 if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
709 }
710
711 auto nb_rows = data.getIndices().size();
712 if (!nb_rows) {
714 }
715 auto nb_base_functions = data.getN().size2();
716
718 CHKERR pythonPtr->blockModes(iD, OP::getCoordsAtGaussPts(),
719 modeCentroids[blockCounter],
720 modeBboxes[blockCounter], blockModes);
721
722 auto nb_integration_pts = getGaussPts().size2();
723 if (blockModes.size2() != 3 * nb_integration_pts) {
724 MOFEM_LOG("WORLD", Sev::error)
725 << "Number of modes does not match number of integration points: "
726 << blockModes.size2() << "!=" << 3 * nb_integration_pts;
727 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "modes/integration points");
728 }
729
730 VectorDouble nf(nb_rows);
731
732 int nb_modes = blockModes.size1();
733 for (auto mode = 0; mode != nb_modes; ++mode) {
734 nf.clear();
735 // get mode
736 auto t_mode = getFTensor1FromPtr<3>(&blockModes(mode, 0));
737 // get element volume
738 const double vol = OP::getMeasure();
739 // get integration weights
740 auto t_w = OP::getFTensor0IntegrationWeight();
741 // get base function gradient on rows
742 auto t_base = data.getFTensor0N();
743 // loop over integration points
744 for (int gg = 0; gg != nb_integration_pts; gg++) {
745
746 // take into account Jacobian
747 const double alpha = t_w * vol;
748 // loop over rows base functions
749 auto t_nf = getFTensor1FromPtr<SPACE_DIM>(nf.data().data());
750 int rr = 0;
751 for (; rr != nb_rows / SPACE_DIM; ++rr) {
752 t_nf(i) += alpha * t_base * t_mode(i);
753 ++t_base;
754 ++t_nf;
755 }
756 for (; rr < nb_base_functions; ++rr)
757 ++t_base;
758 ++t_w; // move to another integration weight
759 ++t_mode; // move to another mode
760 }
761 Vec vec = modeVecs[modeCounter + mode];
762 auto size = data.getIndices().size();
763 auto *indices = data.getIndices().data().data();
764 auto *nf_data = nf.data().data();
765 CHKERR VecSetValues(vec, size, indices, nf_data, ADD_VALUES);
766 }
767
769 }
770
771 private:
772 boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
773 MatrixDouble blockModes;
774 std::vector<std::array<double, 3>> modeCentroids;
775 std::vector<std::array<double, 6>> modeBboxes;
776 int iD;
777 std::vector<SmartPetscObj<Vec>> modeVecs;
778 int blockCounter;
779 int modeCounter;
780 };
781
782 auto solve_bdy = [&]() {
784
785 auto fe_lhs = get_lhs_fe();
786 auto fe_rhs = get_rhs_fe();
787 int block_counter = 0;
788 int mode_counter = 0;
789 for (auto &bc : bcs) {
790 auto id = bc->getMeshsetId();
791 Range ents;
792 CHKERR mField.get_moab().get_entities_by_handle(bc->getMeshset(), ents,
793 true);
794 auto range = boost::make_shared<Range>(ents);
795 auto &pip_rhs = fe_rhs->getOpPtrVector();
796 pip_rhs.push_back(new OpMode("ADJOINT_FIELD", pythonPtr, id, modeVecs,
797 modeCentroids, modeBBoxes, block_counter,
798 mode_counter, range));
799 CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
800 fe_rhs);
801 pip_rhs.pop_back();
802 int nb_modes;
803 CHKERR pythonPtr->numberOfModes(id, nb_modes);
804 ++block_counter;
805 mode_counter += nb_modes;
806 MOFEM_LOG("WORLD", Sev::inform)
807 << "Setting mode block block: " << bc->getName()
808 << " with ID: " << bc->getMeshsetId()
809 << " total modes: " << mode_counter;
810 }
811
812 for (auto &v : modeVecs) {
813 CHKERR VecAssemblyBegin(v);
814 CHKERR VecAssemblyEnd(v);
815 CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
816 CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
817 }
818
819 auto M = createDMMatrix(subset_dm_bdy);
820 fe_lhs->B = M;
821 CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
822 fe_lhs);
823 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
824 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
825
826 auto solver = createKSP(mField.get_comm());
827 CHKERR KSPSetOperators(solver, M, M);
828 CHKERR KSPSetFromOptions(solver);
829 CHKERR KSPSetUp(solver);
830 auto v = createDMVector(subset_dm_bdy);
831 for (auto &f : modeVecs) {
832 CHKERR KSPSolve(solver, f, v);
833 CHKERR VecSwap(f, v);
834 }
835
836 for (auto &v : modeVecs) {
837 CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
838 CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
839 }
840
842 };
843
844 CHKERR solve_bdy();
845
846 auto get_elastic_fe_lhs = [&]() {
847 auto fe = boost::make_shared<DomainEle>(mField);
848 fe->getRuleHook = [](int, int, int p_data) {
849 return 2 * p_data + p_data - 1;
850 };
851 auto &pip = fe->getOpPtrVector();
853 "GEOMETRY");
855 mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
856 return fe;
857 };
858
859 auto get_elastic_fe_rhs = [&]() {
860 auto fe = boost::make_shared<DomainEle>(mField);
861 fe->getRuleHook = [](int, int, int p_data) {
862 return 2 * p_data + p_data - 1;
863 };
864 auto &pip = fe->getOpPtrVector();
866 "GEOMETRY");
868 mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
869 return fe;
870 };
871
872 auto adjoint_gradient_postprocess = [&](auto mode) {
874 auto post_proc_mesh = boost::make_shared<moab::Core>();
875 auto post_proc_begin =
876 boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(mField,
877 post_proc_mesh);
878 auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
879 mField, post_proc_mesh);
880
881 auto geom_vec = boost::make_shared<MatrixDouble>();
882
883 auto post_proc_fe =
884 boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
886 post_proc_fe->getOpPtrVector(), {H1}, "GEOMETRY");
887 post_proc_fe->getOpPtrVector().push_back(
888 new OpCalculateVectorFieldValues<SPACE_DIM>("ADJOINT_FIELD", geom_vec,
889 modeVecs[mode]));
890
892
893 post_proc_fe->getOpPtrVector().push_back(
894
895 new OpPPMap(
896
897 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
898
899 {},
900
901 {{"MODE", geom_vec}},
902
903 {},
904
905 {}
906
907 )
908
909 );
910
912 post_proc_begin->getFEMethod());
913 CHKERR DMoFEMLoopFiniteElements(adjointDM, "ADJOINT_DOMAIN_FE",
914 post_proc_fe);
916 post_proc_begin->getFEMethod());
917
918 CHKERR post_proc_end->writeFile("mode_" + std::to_string(mode) + ".h5m");
919
921 };
922
923 auto solve_domain = [&]() {
925 auto fe_lhs = get_elastic_fe_lhs();
926 auto fe_rhs = get_elastic_fe_rhs();
927 auto v = createDMVector(subset_dm_domain);
928 auto F = vectorDuplicate(v);
929 fe_rhs->f = F;
930
931 auto M = createDMMatrix(subset_dm_domain);
932 fe_lhs->B = M;
933 CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
934 fe_lhs);
935 CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
936 CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
937
938 auto solver = createKSP(mField.get_comm());
939 CHKERR KSPSetOperators(solver, M, M);
940 CHKERR KSPSetFromOptions(solver);
941 CHKERR KSPSetUp(solver);
942
943 int mode_counter = 0;
944 for (auto &f : modeVecs) {
945 CHKERR mField.getInterface<FieldBlas>()->setField(0, "ADJOINT_FIELD");
946 CHKERR DMoFEMMeshToLocalVector(subset_dm_bdy, f, INSERT_VALUES,
947 SCATTER_REVERSE);
948 CHKERR VecZeroEntries(F);
949 CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
950 fe_rhs);
951 CHKERR VecAssemblyBegin(F);
952 CHKERR VecAssemblyEnd(F);
953 CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
954 CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
955 CHKERR KSPSolve(solver, F, v);
956 CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
957 CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
958 CHKERR DMoFEMMeshToLocalVector(subset_dm_domain, v, INSERT_VALUES,
959 SCATTER_REVERSE);
960 auto m = createDMVector(adjointDM);
962 SCATTER_FORWARD);
963 f = m;
964 ++mode_counter;
965 }
967 };
968
969 CHKERR solve_domain();
970
971 for (int i = 0; i < modeVecs.size(); ++i) {
972 CHKERR adjoint_gradient_postprocess(i);
973 }
974
976}
#define FTENSOR_INDEX(DIM, I)
constexpr IntegrationType I
Definition adjoint.cpp:27
constexpr int SPACE_DIM
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
@ F
PetscErrorCode DMMoFEMCreateSubDM(DM subdm, DM dm, const char problem_name[])
Must be called by user to set Sub DM MoFEM data structures.
Definition DMMoFEM.cpp:215
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition DMMoFEM.cpp:238
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition DMMoFEM.cpp:280
auto createDMMatrix(DM dm)
Get smart matrix from DM.
Definition DMMoFEM.hpp:1059
MoFEMErrorCode opFactoryDomainLhs(MoFEM::Interface &m_field, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string field_name, boost::shared_ptr< HookeOps::CommonData > common_ptr, Sev sev)
Assemble domain LHS K factory (LHS first overload) Initializes and pushes operators to assemble the L...
Definition HookeOps.hpp:302
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
Definition Types.hpp:77
FTensor::Tensor1< FTensor::PackPtr< double *, S >, DIM > getFTensor1FromPtr(double *ptr)
Make Tensor1 from pointer.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
OpBaseImpl< PETSC, EdgeEleOp > OpBase
Definition radiation.cpp:29
std::vector< std::array< double, 3 > > modeCentroids
Definition adjoint.cpp:159
std::vector< std::array< double, 6 > > modeBBoxes
Definition adjoint.cpp:160
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
auto save_range

◆ tsSolve()

MoFEMErrorCode Example::tsSolve ( )
private
Examples
plastic.cpp.

Definition at line 834 of file plastic.cpp.

834 {
836
839 ISManager *is_manager = mField.getInterface<ISManager>();
840
841 auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
842
843 auto set_section_monitor = [&](auto solver) {
845 SNES snes;
846 CHKERR TSGetSNES(solver, &snes);
847 CHKERR SNESMonitorSet(snes,
848 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
850 (void *)(snes_ctx_ptr.get()), nullptr);
852 };
853
854 auto create_post_process_elements = [&]() {
855 auto push_vol_ops = [this](auto &pip) {
857 pip, {H1, HDIV}, "GEOMETRY");
858
859 auto [common_plastic_ptr, common_hencky_ptr] =
861 mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
862
863 if (common_hencky_ptr) {
864 if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
865 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
866 }
867
868 return std::make_pair(common_plastic_ptr, common_hencky_ptr);
869 };
870
871 auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
873
874 auto &pip = pp_fe->getOpPtrVector();
875
876 auto [common_plastic_ptr, common_hencky_ptr] = p;
877
879
880 auto x_ptr = boost::make_shared<MatrixDouble>();
881 pip.push_back(
882 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
883 auto u_ptr = boost::make_shared<MatrixDouble>();
884 pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
885
886 if (is_large_strains) {
887
888 pip.push_back(
889
890 new OpPPMap(
891
892 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
893
894 {{"PLASTIC_SURFACE",
895 common_plastic_ptr->getPlasticSurfacePtr()},
896 {"PLASTIC_MULTIPLIER",
897 common_plastic_ptr->getPlasticTauPtr()}},
898
899 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
900
901 {{"GRAD", common_hencky_ptr->matGradPtr},
902 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
903
904 {{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
905 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
906 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
907
908 )
909
910 );
911
912 } else {
913
914 pip.push_back(
915
916 new OpPPMap(
917
918 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
919
920 {{"PLASTIC_SURFACE",
921 common_plastic_ptr->getPlasticSurfacePtr()},
922 {"PLASTIC_MULTIPLIER",
923 common_plastic_ptr->getPlasticTauPtr()}},
924
925 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
926
927 {},
928
929 {{"STRAIN", common_plastic_ptr->mStrainPtr},
930 {"STRESS", common_plastic_ptr->mStressPtr},
931 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
932 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
933
934 )
935
936 );
937 }
938
940 };
941
942 PetscBool post_proc_vol;
943 PetscBool post_proc_skin;
944
945 if constexpr (SPACE_DIM == 2) {
946 post_proc_vol = PETSC_TRUE;
947 post_proc_skin = PETSC_FALSE;
948 } else {
949 post_proc_vol = PETSC_FALSE;
950 post_proc_skin = PETSC_TRUE;
951 }
952 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_vol", &post_proc_vol,
953 PETSC_NULLPTR);
954 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin",
955 &post_proc_skin, PETSC_NULLPTR);
956
957 auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
958 post_proc_vol]() {
959 if (post_proc_vol == PETSC_FALSE)
960 return boost::shared_ptr<PostProcEle>();
961 auto pp_fe = boost::make_shared<PostProcEle>(mField);
963 push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
964 "push_vol_post_proc_ops");
965 return pp_fe;
966 };
967
968 auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
969 post_proc_skin]() {
970 if (post_proc_skin == PETSC_FALSE)
971 return boost::shared_ptr<SkinPostProcEle>();
972
974 auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
975 auto op_side = new OpLoopSide<SideEle>(mField, simple->getDomainFEName(),
976 SPACE_DIM, Sev::verbose);
977 pp_fe->getOpPtrVector().push_back(op_side);
978 CHK_MOAB_THROW(push_vol_post_proc_ops(
979 pp_fe, push_vol_ops(op_side->getOpPtrVector())),
980 "push_vol_post_proc_ops");
981 return pp_fe;
982 };
983
984 return std::make_pair(vol_post_proc(), skin_post_proc());
985 };
986
987 auto scatter_create = [&](auto D, auto coeff) {
989 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
990 ROW, "U", coeff, coeff, is);
991 int loc_size;
992 CHKERR ISGetLocalSize(is, &loc_size);
993 Vec v;
994 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
995 VecScatter scatter;
996 CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
997 return std::make_tuple(SmartPetscObj<Vec>(v),
999 };
1000
1001 boost::shared_ptr<SetPtsData> field_eval_data;
1002 boost::shared_ptr<MatrixDouble> u_field_ptr;
1003
1004 std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
1005 int coords_dim = 3;
1006 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
1007 field_eval_coords.data(), &coords_dim,
1008 &do_eval_field);
1009
1010 boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
1011 scalar_field_ptrs = boost::make_shared<
1012 std::map<std::string, boost::shared_ptr<VectorDouble>>>();
1013 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1014 vector_field_ptrs = boost::make_shared<
1015 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1016 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1017 sym_tensor_field_ptrs = boost::make_shared<
1018 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1019 boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
1020 tensor_field_ptrs = boost::make_shared<
1021 std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
1022
1023 if (do_eval_field) {
1024 auto u_field_ptr = boost::make_shared<MatrixDouble>();
1025 field_eval_data =
1026 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
1027
1028 CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
1029 field_eval_data, simple->getDomainFEName());
1030
1031 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
1032 auto no_rule = [](int, int, int) { return -1; };
1033 auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
1034 field_eval_fe_ptr->getRuleHook = no_rule;
1035
1037 field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, "GEOMETRY");
1038
1039 auto [common_plastic_ptr, common_hencky_ptr] =
1041 mField, "MAT_PLASTIC", field_eval_fe_ptr->getOpPtrVector(), "U",
1042 "EP", "TAU", 1., Sev::inform);
1043
1044 field_eval_fe_ptr->getOpPtrVector().push_back(
1045 new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_field_ptr));
1046
1047 if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
1048 if (is_large_strains) {
1049 scalar_field_ptrs->insert(
1050 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1051 scalar_field_ptrs->insert(
1052 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1053 vector_field_ptrs->insert({"U", u_field_ptr});
1054 sym_tensor_field_ptrs->insert(
1055 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1056 sym_tensor_field_ptrs->insert(
1057 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1058 sym_tensor_field_ptrs->insert(
1059 {"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()});
1060 tensor_field_ptrs->insert({"GRAD", common_hencky_ptr->matGradPtr});
1061 tensor_field_ptrs->insert(
1062 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()});
1063 } else {
1064 scalar_field_ptrs->insert(
1065 {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
1066 scalar_field_ptrs->insert(
1067 {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
1068 vector_field_ptrs->insert({"U", u_field_ptr});
1069 sym_tensor_field_ptrs->insert(
1070 {"STRAIN", common_plastic_ptr->mStrainPtr});
1071 sym_tensor_field_ptrs->insert(
1072 {"STRESS", common_plastic_ptr->mStressPtr});
1073 sym_tensor_field_ptrs->insert(
1074 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
1075 sym_tensor_field_ptrs->insert(
1076 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
1077 }
1078 }
1079 }
1080
1081 auto test_monitor_ptr = boost::make_shared<FEMethod>();
1082
1083 auto set_time_monitor = [&](auto dm, auto solver) {
1085 boost::shared_ptr<Monitor<SPACE_DIM>> monitor_ptr(new Monitor<SPACE_DIM>(
1086 dm, create_post_process_elements(), reactionFe, uXScatter, uYScatter,
1087 uZScatter, field_eval_coords, field_eval_data, scalar_field_ptrs,
1088 vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
1089 boost::shared_ptr<ForcesAndSourcesCore> null;
1090
1091 test_monitor_ptr->postProcessHook = [&]() {
1093
1094 if (atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
1095 test_monitor_ptr->ts_step == 25) {
1096
1097 if (scalar_field_ptrs->at("PLASTIC_MULTIPLIER")->size()) {
1098 auto t_tau =
1099 getFTensor0FromVec(*scalar_field_ptrs->at("PLASTIC_MULTIPLIER"));
1100 MOFEM_LOG("PlasticSync", Sev::inform) << "Eval point tau: " << t_tau;
1101
1102 if (atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
1103 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1104 "atom test %d failed: wrong plastic multiplier value",
1105 atom_test);
1106 }
1107 }
1108
1109 if (vector_field_ptrs->at("U")->size1()) {
1111 auto t_disp =
1112 getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at("U"));
1113 MOFEM_LOG("PlasticSync", Sev::inform) << "Eval point U: " << t_disp;
1114
1115 if (atom_test == 1 && fabs(t_disp(0) - 0.25 / 2.) > 1e-5 ||
1116 fabs(t_disp(1) + 0.0526736) > 1e-5) {
1117 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1118 "atom test %d failed: wrong displacement value",
1119 atom_test);
1120 }
1121 }
1122
1123 if (sym_tensor_field_ptrs->at("PLASTIC_STRAIN")->size1()) {
1124 auto t_plastic_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(
1125 *sym_tensor_field_ptrs->at("PLASTIC_STRAIN"));
1126 MOFEM_LOG("PlasticSync", Sev::inform)
1127 << "Eval point EP: " << t_plastic_strain;
1128
1129 if (atom_test == 1 &&
1130 fabs(t_plastic_strain(0, 0) - 0.221943) > 1e-5 ||
1131 fabs(t_plastic_strain(0, 1)) > 1e-5 ||
1132 fabs(t_plastic_strain(1, 1) + 0.110971) > 1e-5) {
1133 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1134 "atom test %d failed: wrong plastic strain value",
1135 atom_test);
1136 }
1137 }
1138
1139 if (tensor_field_ptrs->at("FIRST_PIOLA")->size1()) {
1140 auto t_piola_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
1141 *tensor_field_ptrs->at("FIRST_PIOLA"));
1142 MOFEM_LOG("PlasticSync", Sev::inform)
1143 << "Eval point Piola stress: " << t_piola_stress;
1144
1145 if (atom_test == 1 && fabs((t_piola_stress(0, 0) - 198.775) /
1146 t_piola_stress(0, 0)) > 1e-5 ||
1147 fabs(t_piola_stress(0, 1)) + fabs(t_piola_stress(1, 0)) +
1148 fabs(t_piola_stress(1, 1)) >
1149 1e-5) {
1150 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
1151 "atom test %d failed: wrong Piola stress value",
1152 atom_test);
1153 }
1154 }
1155 }
1156
1159 };
1160
1161 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
1162 monitor_ptr, null, test_monitor_ptr);
1163
1165 };
1166
1167 auto set_schur_pc = [&](auto solver,
1168 boost::shared_ptr<SetUpSchur> &schur_ptr) {
1170
1171 auto name_prb = simple->getProblemName();
1172
1173 // create sub dm for Schur complement
1174 auto create_schur_dm = [&](SmartPetscObj<DM> base_dm,
1175 SmartPetscObj<DM> &dm_sub) {
1177 dm_sub = createDM(mField.get_comm(), "DMMOFEM");
1178 CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SCHUR");
1179 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
1180 CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
1181 CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
1182 for (auto f : {"U"}) {
1185 }
1186 CHKERR DMSetUp(dm_sub);
1187
1189 };
1190
1191 auto create_block_dm = [&](SmartPetscObj<DM> base_dm,
1192 SmartPetscObj<DM> &dm_sub) {
1194 dm_sub = createDM(mField.get_comm(), "DMMOFEM");
1195 CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "BLOCK");
1196 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
1197 CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
1198 CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
1199#ifdef ADD_CONTACT
1200 for (auto f : {"SIGMA", "EP", "TAU"}) {
1203 }
1204#else
1205 for (auto f : {"EP", "TAU"}) {
1208 }
1209#endif
1210 CHKERR DMSetUp(dm_sub);
1212 };
1213
1214 // Create nested (sub BC) Schur DM
1215 if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
1216
1217 SmartPetscObj<DM> dm_schur;
1218 CHKERR create_schur_dm(simple->getDM(), dm_schur);
1219 SmartPetscObj<DM> dm_block;
1220 CHKERR create_block_dm(simple->getDM(), dm_block);
1221
1222#ifdef ADD_CONTACT
1223
1224 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1225 auto block_mat_data = createBlockMatStructure(
1226 simple->getDM(),
1227
1228 {
1229
1230 {simple->getDomainFEName(),
1231
1232 {{"U", "U"},
1233 {"SIGMA", "SIGMA"},
1234 {"U", "SIGMA"},
1235 {"SIGMA", "U"},
1236 {"EP", "EP"},
1237 {"TAU", "TAU"},
1238 {"U", "EP"},
1239 {"EP", "U"},
1240 {"EP", "TAU"},
1241 {"TAU", "EP"},
1242 {"TAU", "U"}
1243
1244 }},
1245
1246 {simple->getBoundaryFEName(),
1247
1248 {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
1249
1250 }}
1251
1252 }
1253
1254 );
1255
1257
1258 {dm_schur, dm_block}, block_mat_data,
1259
1260 {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
1261
1262 );
1263 };
1264
1265#else
1266
1267 auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
1268 auto block_mat_data =
1270
1271 {{simple->getDomainFEName(),
1272
1273 {{"U", "U"},
1274 {"EP", "EP"},
1275 {"TAU", "TAU"},
1276 {"U", "EP"},
1277 {"EP", "U"},
1278 {"EP", "TAU"},
1279 {"TAU", "U"},
1280 {"TAU", "EP"}
1281
1282 }}}
1283
1284 );
1285
1287
1288 {dm_schur, dm_block}, block_mat_data,
1289
1290 {"EP", "TAU"}, {nullptr, nullptr}, false
1291
1292 );
1293 };
1294
1295#endif
1296
1297 auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
1298 CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
1299
1300 auto block_is = getDMSubData(dm_block)->getSmartRowIs();
1301 auto ao_schur = getDMSubData(dm_schur)->getSmartRowMap();
1302
1303 // Indices has to be map fro very to level, while assembling Schur
1304 // complement.
1305 schur_ptr =
1306 SetUpSchur::createSetUpSchur(mField, dm_schur, block_is, ao_schur);
1307 CHKERR schur_ptr->setUp(solver);
1308 }
1309
1311 };
1312
1313 auto dm = simple->getDM();
1314 auto D = createDMVector(dm);
1315 auto DD = vectorDuplicate(D);
1316 uXScatter = scatter_create(D, 0);
1317 uYScatter = scatter_create(D, 1);
1318 if constexpr (SPACE_DIM == 3)
1319 uZScatter = scatter_create(D, 2);
1320
1321 auto create_solver = [pip_mng]() {
1322 if (is_quasi_static == PETSC_TRUE)
1323 return pip_mng->createTSIM();
1324 else
1325 return pip_mng->createTSIM2();
1326 };
1327
1328 auto solver = create_solver();
1329
1330 auto active_pre_lhs = []() {
1332 std::fill(PlasticOps::CommonData::activityData.begin(),
1335 };
1336
1337 auto active_post_lhs = [&]() {
1339 auto get_iter = [&]() {
1340 SNES snes;
1341 CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
1342 int iter;
1343 CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
1344 "Can not get iter");
1345 return iter;
1346 };
1347
1348 auto iter = get_iter();
1349 if (iter >= 0) {
1350
1351 std::array<int, 5> activity_data;
1352 std::fill(activity_data.begin(), activity_data.end(), 0);
1353 MPI_Allreduce(PlasticOps::CommonData::activityData.data(),
1354 activity_data.data(), activity_data.size(), MPI_INT,
1355 MPI_SUM, mField.get_comm());
1356
1357 int &active_points = activity_data[0];
1358 int &avtive_full_elems = activity_data[1];
1359 int &avtive_elems = activity_data[2];
1360 int &nb_points = activity_data[3];
1361 int &nb_elements = activity_data[4];
1362
1363 if (nb_points) {
1364
1365 double proc_nb_points =
1366 100 * static_cast<double>(active_points) / nb_points;
1367 double proc_nb_active =
1368 100 * static_cast<double>(avtive_elems) / nb_elements;
1369 double proc_nb_full_active = 100;
1370 if (avtive_elems)
1371 proc_nb_full_active =
1372 100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1373
1374 MOFEM_LOG_C("PLASTICITY", Sev::inform,
1375 "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
1376 "elements %d "
1377 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1378 iter, nb_points, active_points, proc_nb_points,
1379 avtive_elems, proc_nb_active, avtive_full_elems,
1380 proc_nb_full_active, iter);
1381 }
1382 }
1383
1385 };
1386
1387 auto add_active_dofs_elem = [&](auto dm) {
1389 auto fe_pre_proc = boost::make_shared<FEMethod>();
1390 fe_pre_proc->preProcessHook = active_pre_lhs;
1391 auto fe_post_proc = boost::make_shared<FEMethod>();
1392 fe_post_proc->postProcessHook = active_post_lhs;
1393 auto ts_ctx_ptr = getDMTsCtx(dm);
1394 ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
1395 ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
1397 };
1398
1399 auto set_essential_bc = [&](auto dm, auto solver) {
1401 // This is low level pushing finite elements (pipelines) to solver
1402
1403 auto pre_proc_ptr = boost::make_shared<FEMethod>();
1404 auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
1405 auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
1406
1407 // Add boundary condition scaling
1408 auto disp_time_scale = boost::make_shared<TimeScale>();
1409
1410 auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
1412 {disp_time_scale}, false);
1413 return hook;
1414 };
1415 pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
1416
1417 auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
1420 mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
1422 mField, post_proc_rhs_ptr, 1.)();
1424 };
1425 auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
1427 mField, post_proc_lhs_ptr, 1.);
1428 };
1429 post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
1430
1431 auto ts_ctx_ptr = getDMTsCtx(dm);
1432 ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
1433 ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
1434 ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
1435 post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
1436 ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
1437
1439 };
1440
1441 auto B = createDMMatrix(dm);
1442 if (is_quasi_static == PETSC_FALSE) {
1443 CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
1444 } else {
1445 CHKERR TSSetI2Jacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
1446 }
1447 if (is_quasi_static == PETSC_TRUE) {
1448 CHKERR TSSetSolution(solver, D);
1449 } else {
1450 CHKERR TS2SetSolution(solver, D, DD);
1451 }
1452 CHKERR set_section_monitor(solver);
1453 CHKERR set_time_monitor(dm, solver);
1454 CHKERR TSSetFromOptions(solver);
1455
1456 CHKERR add_active_dofs_elem(dm);
1457 boost::shared_ptr<SetUpSchur> schur_ptr;
1458 CHKERR set_schur_pc(solver, schur_ptr);
1459 CHKERR set_essential_bc(dm, solver);
1460
1461 MOFEM_LOG_CHANNEL("TIMER");
1462 MOFEM_LOG_TAG("TIMER", "timer");
1463 if (set_timer)
1464 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1465 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
1466 CHKERR TSSetUp(solver);
1467 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
1468 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
1469 CHKERR TSSolve(solver, NULL);
1470 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
1471
1473}
#define CHK_MOAB_THROW(err, msg)
Check error code of MoAB function and throw MoFEM exception.
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition DMMoFEM.hpp:1144
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *ctx)
Sens monitor printing residual field by field.
Definition SnesCtx.cpp:592
auto getDMSubData(DM dm)
Get sub problem data structure.
Definition DMMoFEM.hpp:1160
FTensor::Tensor2< FTensor::PackPtr< double *, 1 >, Tensor_Dim1, Tensor_Dim2 > getFTensor2FromMat(MatrixDouble &data)
Get tensor rank 2 (matrix) form data matrix.
boost::shared_ptr< BlockStructure > createBlockMatStructure(DM dm, SchurFEOpsFEandFields schur_fe_op_vec)
Create a Mat Diag Blocks object.
Definition Schur.cpp:1082
boost::shared_ptr< NestSchurData > createSchurNestedMatrixStruture(std::pair< SmartPetscObj< DM >, SmartPetscObj< DM > > dms, boost::shared_ptr< BlockStructure > block_mat_data_ptr, std::vector< std::string > fields_names, std::vector< boost::shared_ptr< Range > > field_ents, bool add_preconditioner_block)
Get the Schur Nest Mat Array object.
Definition Schur.cpp:2343
static auto getFTensor2SymmetricFromMat(ublas::matrix< T, L, A > &data)
Get symmetric tensor rank 2 (matrix) form data matrix.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode DMMoFEMSetNestSchurData(DM dm, boost::shared_ptr< NestSchurData >)
Definition DMMoFEM.cpp:1554
auto getDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition DMMoFEM.hpp:1130
auto createCommonPlasticOps(MoFEM::Interface &m_field, std::string block_name, boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pip, std::string u, std::string ep, std::string tau, double scale, Sev sev)
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition plastic.cpp:239
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition plastic.cpp:237
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition plastic.cpp:238
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
static std::array< int, 5 > activityData
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field)

Friends And Related Symbol Documentation

◆ TSPrePostProc

friend struct TSPrePostProc
friend

Definition at line 434 of file dynamic_first_order_con_law.cpp.

Member Data Documentation

◆ adjointDM

SmartPetscObj<DM> Example::adjointDM
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 156 of file adjoint.cpp.

◆ approxFunction

ApproxFieldFunction< FIELD_DIM > Example::approxFunction
staticprivate
Initial value:

Definition at line 61 of file approximaton.cpp.

◆ approxGradVals

boost::shared_ptr<MatrixDouble> Example::approxGradVals
private

Definition at line 63 of file radiation.cpp.

◆ approxVals

boost::shared_ptr<VectorDouble> Example::approxVals
private

Definition at line 62 of file radiation.cpp.

◆ aveMaxMin

std::array< double, Example::BoundingBox::LAST_BB > Example::aveMaxMin
staticprivate
Examples
phase.cpp.

Definition at line 110 of file phase.cpp.

◆ base

◆ boundaryMarker

boost::shared_ptr< std::vector< unsigned char > > Example::boundaryMarker
private
Examples
helmholtz.cpp, and phase.cpp.

Definition at line 42 of file helmholtz.cpp.

◆ commonDataPtr

boost::shared_ptr< CommonData > Example::commonDataPtr
private

Definition at line 42 of file integration.cpp.

◆ domianLhsFEPtr

boost::shared_ptr<FEMethod> Example::domianLhsFEPtr
private
Examples
shallow_wave.cpp.

Definition at line 355 of file shallow_wave.cpp.

◆ domianRhsFEPtr

boost::shared_ptr<FEMethod> Example::domianRhsFEPtr
private
Examples
shallow_wave.cpp.

Definition at line 356 of file shallow_wave.cpp.

◆ ePS

SmartPetscObj<EPS> Example::ePS
private
Examples
eigen_elastic.cpp.

Definition at line 69 of file eigen_elastic.cpp.

◆ fieldOrder

int Example::fieldOrder = 2
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 153 of file adjoint.cpp.

◆ focalIndex

int Example::focalIndex
staticprivate
Examples
phase.cpp.

Definition at line 112 of file phase.cpp.

◆ iI

std::vector< MatrixInt > Example::iI
staticprivate
Examples
phase.cpp.

Definition at line 109 of file phase.cpp.

◆ initialGeometry

SmartPetscObj<Vec> Example::initialGeometry
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 161 of file adjoint.cpp.

◆ K

SmartPetscObj<Mat> Example::K
private
Examples
eigen_elastic.cpp.

Definition at line 68 of file eigen_elastic.cpp.

◆ kspElastic

SmartPetscObj<KSP> Example::kspElastic
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 155 of file adjoint.cpp.

◆ M

SmartPetscObj<Mat> Example::M
private
Examples
eigen_elastic.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 67 of file eigen_elastic.cpp.

◆ matDPtr

boost::shared_ptr<MatrixDouble> Example::matDPtr
private
Examples
eigen_elastic.cpp.

Definition at line 65 of file eigen_elastic.cpp.

◆ meshVolumeAndCount

std::array<double, 2> Example::meshVolumeAndCount = {0, 0}
inlinestatic
Examples
plastic.cpp.

Definition at line 223 of file plastic.cpp.

223{0, 0};

◆ mField

◆ modeBBoxes

std::vector<std::array<double, 6> > Example::modeBBoxes
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 160 of file adjoint.cpp.

◆ modeCentroids

std::vector<std::array<double, 3> > Example::modeCentroids
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 159 of file adjoint.cpp.

◆ modeVecs

std::vector<SmartPetscObj<Vec> > Example::modeVecs
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 158 of file adjoint.cpp.

◆ pinchNodes

Range Example::pinchNodes
private
Examples
heat_method.cpp.

Definition at line 54 of file heat_method.cpp.

◆ pythonPtr

boost::shared_ptr<ObjectiveFunctionData> Example::pythonPtr
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 157 of file adjoint.cpp.

◆ reactionFe

boost::shared_ptr<DomainEle> Example::reactionFe
private
Examples
plastic.cpp.

Definition at line 235 of file plastic.cpp.

◆ rigidBodyMotion

std::array<SmartPetscObj<Vec>, 6> Example::rigidBodyMotion
private
Examples
eigen_elastic.cpp.

Definition at line 71 of file eigen_elastic.cpp.

◆ rZ

std::vector< double > Example::rZ
staticprivate
Examples
phase.cpp.

Definition at line 108 of file phase.cpp.

◆ savitzkyGolayNormalisation

int Example::savitzkyGolayNormalisation
staticprivate
Examples
phase.cpp.

Definition at line 113 of file phase.cpp.

◆ savitzkyGolayWeights

const int * Example::savitzkyGolayWeights
staticprivate
Examples
phase.cpp.

Definition at line 114 of file phase.cpp.

◆ simpleInterface

Simple * Example::simpleInterface
private
Examples
heat_method.cpp, helmholtz.cpp, phase.cpp, and plot_base.cpp.

Definition at line 41 of file helmholtz.cpp.

◆ space

FieldSpace Example::space
private
Examples
plot_base.cpp.

Definition at line 69 of file plot_base.cpp.

◆ uXScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Example::uXScatter
private
Examples
plastic.cpp.

Definition at line 237 of file plastic.cpp.

◆ uYScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Example::uYScatter
private
Examples
plastic.cpp.

Definition at line 238 of file plastic.cpp.

◆ uZScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Example::uZScatter
private
Examples
plastic.cpp.

Definition at line 239 of file plastic.cpp.

◆ vectorFieldPtr

boost::shared_ptr< MatrixDouble > Example::vectorFieldPtr = nullptr
private
Examples
mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 87 of file elastic.cpp.


The documentation for this struct was generated from the following files: