v0.13.2
Loading...
Searching...
No Matches
Classes | Public Member Functions | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes | Static Private Attributes | List of all members
Example Struct Reference

[Example] More...

Collaboration diagram for Example:
[legend]

Classes

struct  BoundaryOp
 
struct  CommonData
 [Example] More...
 
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  PlasticityTimeScale
 

Public Member Functions

 Example (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProblem ()
 [Run problem] More...
 
 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 ()
 

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] More...
 
MoFEMErrorCode createCommonData ()
 [Set up problem] More...
 
MoFEMErrorCode bC ()
 [Create common data] More...
 
MoFEMErrorCode OPs ()
 [Boundary condition] More...
 
MoFEMErrorCode tsSolve ()
 
MoFEMErrorCode readMesh ()
 [run problem] More...
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode boundaryCondition ()
 [Set up problem] More...
 
MoFEMErrorCode assembleSystem ()
 [Applying essential BC] More...
 
MoFEMErrorCode solveSystem ()
 [Push operators to pipeline] More...
 
MoFEMErrorCode outputResults ()
 [Solve] More...
 
MoFEMErrorCode checkResults ()
 [Postprocess results] More...
 
MoFEMErrorCode setUp ()
 [Run all] More...
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode setFieldValues ()
 [Create common data] More...
 
MoFEMErrorCode pushOperators ()
 [Set density distribution] More...
 
MoFEMErrorCode integrateElements ()
 [Push operators to pipeline] More...
 
MoFEMErrorCode postProcess ()
 [Integrate] More...
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 [Set up problem] More...
 
MoFEMErrorCode createCommonData ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 
MoFEMErrorCode checkResults ()
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode assembleSystemIntensity ()
 [Calculate flux on boundary] More...
 
MoFEMErrorCode assembleSystemFlux ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode calculateFlux (double &calc_flux)
 [Set up problem] More...
 
MoFEMErrorCode outputResults (const int i)
 [Solve] More...
 
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] More...
 
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 addMatBlockOps (boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
 
MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode createCommonData ()
 
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 createCommonData ()
 
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 ()
 

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< PlasticOps::CommonDatacommonPlasticDataPtr
 
boost::shared_ptr< HenckyOps::CommonDatacommonHenckyDataPtr
 
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
 
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
 
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
 
SimplesimpleInterface
 
boost::shared_ptr< CommonDatacommonDataPtr
 
FieldApproximationBase base
 
FieldSpace space
 
boost::shared_ptr< VectorDouble > approxVals
 
boost::shared_ptr< MatrixDouble > approxGradVals
 
Range pinchNodes
 
boost::shared_ptr< MatrixDouble > matDPtr
 
SmartPetscObj< Mat > M
 
SmartPetscObj< Mat > K
 
SmartPetscObj< EPS > ePS
 
std::array< SmartPetscObj< Vec >, 6 > rigidBodyMotion
 
boost::shared_ptr< MatrixDouble > matGradPtr
 
boost::shared_ptr< MatrixDouble > matStrainPtr
 
boost::shared_ptr< MatrixDouble > matStressPtr
 
boost::shared_ptr< MatrixDouble > matAccelerationPtr
 
boost::shared_ptr< MatrixDouble > matInertiaPtr
 
boost::shared_ptr< MatrixDouble > matTangentPtr
 
boost::shared_ptr< FEMethoddomianLhsFEPtr
 
boost::shared_ptr< FEMethoddomianRhsFEPtr
 

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
 

Detailed Description

[Example]

Examples
eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, phase.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 141 of file plastic.cpp.

Member Enumeration Documentation

◆ 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/13]

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

Definition at line 143 of file plastic.cpp.

143: mField(m_field) {}
MoFEM::Interface & mField
Definition: plastic.cpp:148

◆ Example() [2/13]

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

Definition at line 39 of file helmholtz.cpp.

39: mField(m_field) {}

◆ Example() [3/13]

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

Definition at line 25 of file integration.cpp.

25: mField(m_field) {}

◆ Example() [4/13]

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

Definition at line 50 of file plot_base.cpp.

50: mField(m_field) {}

◆ Example() [5/13]

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

Definition at line 81 of file phase.cpp.

81: mField(m_field) {}

◆ Example() [6/13]

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

Definition at line 53 of file approximaton.cpp.

53: mField(m_field) {}

◆ Example() [7/13]

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

Definition at line 45 of file radiation.cpp.

45: mField(m_field) {}

◆ Example() [8/13]

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

Definition at line 46 of file heat_method.cpp.

46: mField(m_field) {}

◆ Example() [9/13]

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

Definition at line 64 of file elastic.cpp.

64: mField(m_field) {}

◆ Example() [10/13]

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

Definition at line 49 of file eigen_elastic.cpp.

49: mField(m_field) {}

◆ Example() [11/13]

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

Definition at line 60 of file nonlinear_elastic.cpp.

60: mField(m_field) {}

◆ Example() [12/13]

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

Definition at line 58 of file nonlinear_dynamic_elastic.cpp.

58: mField(m_field) {}

◆ Example() [13/13]

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

Definition at line 339 of file shallow_wave.cpp.

339: mField(m_field) {}

Member Function Documentation

◆ addMatBlockOps()

MoFEMErrorCode Example::addMatBlockOps ( boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &  pipeline,
std::string  field_name,
std::string  block_name,
boost::shared_ptr< MatrixDouble >  mat_D_Ptr,
Sev  sev 
)
private

[Calculate elasticity tensor]

[Calculate elasticity tensor]

Definition at line 86 of file elastic.cpp.

89 {
91
92 struct OpMatBlocks : public DomainEleOp {
93 OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
94 double bulk_modulus_K, double shear_modulus_G,
95 MoFEM::Interface &m_field, Sev sev,
96 std::vector<const CubitMeshSets *> meshset_vec_ptr)
97 : DomainEleOp(field_name, DomainEleOp::OPROW), matDPtr(m),
98 bulkModulusKDefault(bulk_modulus_K),
99 shearModulusGDefault(shear_modulus_G) {
100 std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
101 CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
102 "Can not get data from block");
103 }
104
105 MoFEMErrorCode doWork(int side, EntityType type,
108
109 for (auto &b : blockData) {
110
111 if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
112 CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
114 }
115 }
116
117 CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
119 }
120
121 private:
122 boost::shared_ptr<MatrixDouble> matDPtr;
123
124 struct BlockData {
125 double bulkModulusK;
126 double shearModulusG;
127 Range blockEnts;
128 };
129
130 double bulkModulusKDefault;
131 double shearModulusGDefault;
132 std::vector<BlockData> blockData;
133
135 extractBlockData(MoFEM::Interface &m_field,
136 std::vector<const CubitMeshSets *> meshset_vec_ptr,
137 Sev sev) {
139
140 for (auto m : meshset_vec_ptr) {
141 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
142 std::vector<double> block_data;
143 CHKERR m->getAttributes(block_data);
144 if (block_data.size() != 2) {
145 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
146 "Expected that block has two attribute");
147 }
148 auto get_block_ents = [&]() {
149 Range ents;
150 CHKERR
151 m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
152 return ents;
153 };
154
155 double young_modulus = block_data[0];
156 double poisson_ratio = block_data[1];
157 double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
158 double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
159
160 MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
161 << "E = " << young_modulus << " nu = " << poisson_ratio;
162
163 blockData.push_back(
164 {bulk_modulus_K, shear_modulus_G, get_block_ents()});
165 }
166 MOFEM_LOG_CHANNEL("WORLD");
168 }
169
170 MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
171 double bulk_modulus_K, double shear_modulus_G) {
173 //! [Calculate elasticity tensor]
174 auto set_material_stiffness = [&]() {
180 double A = (SPACE_DIM == 2)
181 ? 2 * shear_modulus_G /
182 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
183 : 1;
184 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
185 t_D(i, j, k, l) =
186 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
187 A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
188 t_kd(k, l);
189 };
190 //! [Calculate elasticity tensor]
191 constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
192 mat_D_ptr->resize(size_symm * size_symm, 1);
193 set_material_stiffness();
195 }
196 };
197
198 pipeline.push_back(new OpMatBlocks(
200
201 // Get blockset using regular expression
203
204 (boost::format("%s(.*)") % block_name).str()
205
206 ))
207
208 ));
209
211}
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition: LogManager.hpp:359
constexpr int SPACE_DIM
Kronecker Delta class symmetric.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:595
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
constexpr double shear_modulus_G
Definition: elastic.cpp:60
constexpr double bulk_modulus_K
Definition: elastic.cpp:59
FTensor::Index< 'm', SPACE_DIM > m
constexpr auto t_kd
SeverityLevel
Severity levels.
Definition: LogManager.hpp:33
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
MoFEMErrorCode getCubitMeshsetPtr(const int ms_id, const CubitBCType cubit_bc_type, const CubitMeshSets **cubit_meshset_ptr) const
get cubit meshset
FTensor::Index< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
constexpr AssemblyType A
double young_modulus
Definition: plastic.cpp:99
constexpr auto size_symm
Definition: plastic.cpp:33
constexpr auto field_name
virtual moab::Interface & get_moab()=0
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ assembleSystem() [1/9]

MoFEMErrorCode Example::assembleSystem ( )
private

[Applying essential BC]

[Boundary condition]

[Push operators to pipeline]

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

Definition at line 154 of file helmholtz.cpp.

154 {
157
158 double k = 90;
159 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-k", &k, PETSC_NULL);
160
161 auto beta = [](const double, const double, const double) { return -1; };
162 auto k2 = [k](const double, const double, const double) { return pow(k, 2); };
163 auto kp = [k](const double, const double, const double) { return k; };
164 auto km = [k](const double, const double, const double) { return -k; };
165 auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
166
167 auto set_domain = [&]() {
169 auto det_ptr = boost::make_shared<VectorDouble>();
170 auto jac_ptr = boost::make_shared<MatrixDouble>();
171 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
172 pipeline_mng->getOpDomainLhsPipeline().push_back(
173 new OpCalculateHOJac<2>(jac_ptr));
174 pipeline_mng->getOpDomainLhsPipeline().push_back(
175 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
176 pipeline_mng->getOpDomainLhsPipeline().push_back(
177 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
178 pipeline_mng->getOpDomainLhsPipeline().push_back(
180
181 pipeline_mng->getOpDomainLhsPipeline().push_back(
182 new OpSetBc("P_REAL", true, boundaryMarker));
183
184 pipeline_mng->getOpDomainLhsPipeline().push_back(
185 new OpDomainGradGrad("P_REAL", "P_REAL", beta));
186 pipeline_mng->getOpDomainLhsPipeline().push_back(
187 new OpDomainGradGrad("P_IMAG", "P_IMAG", beta));
188
189 pipeline_mng->getOpDomainLhsPipeline().push_back(
190 new OpDomainMass("P_REAL", "P_REAL", k2));
191 pipeline_mng->getOpDomainLhsPipeline().push_back(
192 new OpDomainMass("P_IMAG", "P_IMAG", k2));
193
194 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
195
198 };
199
200 auto set_boundary = [&]() {
202 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
203 new OpSetBc("P_REAL", true, boundaryMarker));
204 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
205 new OpBoundaryMass("P_IMAG", "P_REAL", kp));
206 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
207 new OpBoundaryMass("P_REAL", "P_IMAG", km));
208 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
209
210 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
211 new OpSetBc("P_REAL", false, boundaryMarker));
212 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
213 new OpBoundaryMass("P_REAL", "P_REAL", beta));
214 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
215
216 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
217 new OpSetBc("P_REAL", false, boundaryMarker));
218 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
219 new OpBoundarySource("P_REAL", beta));
220 pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpUnSetBc("P_REAL"));
221
225 };
226
227 CHKERR set_domain();
228 CHKERR set_boundary();
229
231}
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
@ H1
continuous field
Definition: definitions.h:85
auto integration_rule
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:27
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:33
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:73
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:164
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)

◆ assembleSystem() [2/9]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [3/9]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [4/9]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [5/9]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [6/9]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [7/9]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [8/9]

MoFEMErrorCode Example::assembleSystem ( )
private

◆ assembleSystem() [9/9]

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
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
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.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]

◆ bC() [1/2]

MoFEMErrorCode Example::bC ( )
private

[Create common data]

[Boundary condition]

Examples
plastic.cpp.

Definition at line 397 of file plastic.cpp.

397 {
399
401 auto bc_mng = mField.getInterface<BcManager>();
402 auto prb_mng = mField.getInterface<ProblemsManager>();
403
404 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
405 "U", 0, 0);
406 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
407 "U", 1, 1);
408 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
409 "U", 2, 2);
410 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
411 "REMOVE_ALL", "U", 0, 3);
412
413 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
414 simple->getProblemName(), "U");
415
416 auto &bc_map = bc_mng->getBcMapByBlockName();
417 boundaryMarker = bc_mng->getMergedBlocksMarker(vector<string>{"FIX_"});
418
419 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "REACTION",
420 "U", 0, 3);
421
422 for (auto bc : bc_map)
423 MOFEM_LOG("EXAMPLE", Sev::verbose) << "Marker " << bc.first;
424
425 // OK. We have problem with GMesh, it adding empty characters at the end of
426 // block. So first block is search by regexp. popMarkDOFsOnEntities should
427 // work with regexp.
428 std::string reaction_block_set;
429 for (auto bc : bc_map) {
430 if (bc_mng->checkBlock(bc, "REACTION")) {
431 reaction_block_set = bc.first;
432 break;
433 }
434 }
435
436 if (auto bc = bc_mng->popMarkDOFsOnEntities(reaction_block_set)) {
437 reactionMarker = bc->getBcMarkersPtr();
438
439 // Only take reaction from nodes
440 Range nodes;
441 CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes, true);
442 CHKERR prb_mng->markDofs(simple->getProblemName(), ROW,
443 ProblemsManager::MarkOP::AND, nodes,
445
446 } else {
447 MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
448 }
449
450 if (!reactionMarker) {
451 MOFEM_LOG("EXAMPLE", Sev::warning) << "REACTION blockset does not exist";
452 }
453
455}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
@ ROW
Definition: definitions.h:123
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
boost::shared_ptr< std::vector< unsigned char > > reactionMarker
Definition: plastic.cpp:165
Simple interface for fast problem set-up.
Definition: BcManager.hpp:24
Definition of the displacement bc data structure.
Definition: BCData.hpp:72
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:27

◆ bC() [2/2]

MoFEMErrorCode Example::bC ( )
private

◆ boundaryCondition() [1/9]

MoFEMErrorCode Example::boundaryCondition ( )
private

[Set up problem]

[Create common data]

[Applying essential BC]

[Boundary condition]

[Define gravity vector]

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

Definition at line 106 of file helmholtz.cpp.

106 {
108
109 auto get_ents_on_mesh_skin = [&]() {
110 Range boundary_entities;
112 std::string entity_name = it->getName();
113 if (entity_name.compare(0, 2, "BC") == 0) {
114 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
115 boundary_entities, true);
116 }
117 }
118 // Add vertices to boundary entities
119 Range boundary_vertices;
120 CHKERR mField.get_moab().get_connectivity(boundary_entities,
121 boundary_vertices, true);
122 boundary_entities.merge(boundary_vertices);
123
124 return boundary_entities;
125 };
126
127 auto mark_boundary_dofs = [&](Range &&skin_edges) {
128 auto problem_manager = mField.getInterface<ProblemsManager>();
129 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
130 problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
131 skin_edges, *marker_ptr);
132 return marker_ptr;
133 };
134
135 auto remove_dofs_from_problem = [&](Range &&skin_edges) {
137 auto problem_manager = mField.getInterface<ProblemsManager>();
138 CHKERR problem_manager->removeDofsOnEntities(
139 simpleInterface->getProblemName(), "P_IMAG", skin_edges, 0, 1);
141 };
142
143 // Get global local vector of marked DOFs. Is global, since is set for all
144 // DOFs on processor. Is local since only DOFs on processor are in the
145 // vector. To access DOFs use local indices.
146 boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
147 CHKERR remove_dofs_from_problem(get_ents_on_mesh_skin());
148
150}
@ BLOCKSET
Definition: definitions.h:148
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
Simple * simpleInterface
Definition: helmholtz.cpp:45
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:348

◆ boundaryCondition() [2/9]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [3/9]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [4/9]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [5/9]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [6/9]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [7/9]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [8/9]

MoFEMErrorCode Example::boundaryCondition ( )
private

◆ boundaryCondition() [9/9]

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.

◆ checkResults() [1/11]

MoFEMErrorCode Example::checkResults ( )
private

[Postprocess results]

[Postprocessing results]

[Print results]

[Check results]

[Test example]

[Check]

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

Definition at line 299 of file helmholtz.cpp.

299 {
302
303 auto dm = simpleInterface->getDM();
304 auto D = smartCreateDMVector(dm);
305 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
306 double nrm2;
307 CHKERR VecNorm(D, NORM_2, &nrm2);
308 MOFEM_LOG("WORLD", Sev::inform)
309 << std::setprecision(9) << "Solution norm " << nrm2;
310
311 PetscBool test_flg = PETSC_FALSE;
312 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
313 if (test_flg) {
314 constexpr double regression_test = 97.261672;
315 constexpr double eps = 1e-6;
316 if (std::abs(nrm2 - regression_test) / regression_test > eps)
317 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
318 "Not converged solution");
319 }
321}
static const double eps
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:511
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:995
double D
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:671

◆ 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/10]

MoFEMErrorCode Example::createCommonData ( )
private

[Set up problem]

[Set integration rule]

[Create common data]

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

Definition at line 261 of file plastic.cpp.

261 {
263
264 auto get_command_line_parameters = [&]() {
266 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
267 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
268 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
269 &young_modulus, PETSC_NULL);
270 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
271 &poisson_ratio, PETSC_NULL);
272 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
273 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
274 PETSC_NULL);
275 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
276 PETSC_NULL);
277 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn0", &cn0, PETSC_NULL);
278 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn1", &cn1, PETSC_NULL);
279 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-zeta", &zeta, PETSC_NULL);
280 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
281 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
282 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
283 &is_large_strains, PETSC_NULL);
284 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-set_timer", &set_timer,
285 PETSC_NULL);
286
287 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
288 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
289 PETSC_NULL);
290
291 MOFEM_LOG("EXAMPLE", Sev::inform) << "Young modulus " << young_modulus;
292 MOFEM_LOG("EXAMPLE", Sev::inform) << "Poisson ratio " << poisson_ratio;
293 MOFEM_LOG("EXAMPLE", Sev::inform) << "Yield stress " << sigmaY;
294 MOFEM_LOG("EXAMPLE", Sev::inform) << "Hardening " << H;
295 MOFEM_LOG("EXAMPLE", Sev::inform) << "Viscous hardening " << visH;
296 MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation yield stress " << Qinf;
297 MOFEM_LOG("EXAMPLE", Sev::inform) << "Saturation exponent " << b_iso;
298 MOFEM_LOG("EXAMPLE", Sev::inform) << "cn0 " << cn0;
299 MOFEM_LOG("EXAMPLE", Sev::inform) << "cn1 " << cn1;
300 MOFEM_LOG("EXAMPLE", Sev::inform) << "zeta " << zeta;
301
302 MOFEM_LOG("EXAMPLE", Sev::inform) << "order " << order;
303 MOFEM_LOG("EXAMPLE", Sev::inform) << "geom order " << geom_order;
304
305 PetscBool is_scale = PETSC_TRUE;
306 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_scale", &is_scale,
307 PETSC_NULL);
308 if (is_scale) {
311 rho *= scale;
312 sigmaY *= scale;
313 H *= scale;
314 Qinf *= scale;
315 visH *= scale;
316
317 MOFEM_LOG("EXAMPLE", Sev::inform)
318 << "Scaled Young modulus " << young_modulus;
319 MOFEM_LOG("EXAMPLE", Sev::inform)
320 << "Scaled Poisson ratio " << poisson_ratio;
321 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Yield stress " << sigmaY;
322 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Hardening " << H;
323 MOFEM_LOG("EXAMPLE", Sev::inform) << "Scaled Viscous hardening " << visH;
324 MOFEM_LOG("EXAMPLE", Sev::inform)
325 << "Scaled Saturation yield stress " << Qinf;
326 }
327
329 };
330
331 auto set_matrial_stiffness = [&]() {
338 const double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
339 const double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
340
341 // Plane stress or when 1, plane strain or 3d
342 const double A = (SPACE_DIM == 2)
343 ? 2 * shear_modulus_G /
344 (bulk_modulus_K + (4. / 3.) * shear_modulus_G)
345 : 1;
346
347 auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
348 *commonPlasticDataPtr->mDPtr);
349 auto t_D_axiator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
350 *commonPlasticDataPtr->mDPtr_Axiator);
351 auto t_D_deviator = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(
352 *commonPlasticDataPtr->mDPtr_Deviator);
353
354 constexpr double third = boost::math::constants::third<double>();
355 t_D_axiator(i, j, k, l) = A *
356 (bulk_modulus_K - (2. / 3.) * shear_modulus_G) *
357 t_kd(i, j) * t_kd(k, l);
358 t_D_deviator(i, j, k, l) =
359 2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.);
360 t_D(i, j, k, l) = t_D_axiator(i, j, k, l) + t_D_deviator(i, j, k, l);
361
363 };
364
365 auto make_d_mat = []() {
366 return boost::make_shared<MatrixDouble>(size_symm * size_symm, 1);
367 };
368
369 commonPlasticDataPtr = boost::make_shared<PlasticOps::CommonData>();
370 commonPlasticDataPtr->mDPtr = make_d_mat();
371 commonPlasticDataPtr->mDPtr_Axiator = make_d_mat();
372 commonPlasticDataPtr->mDPtr_Deviator = make_d_mat();
373
374 commonPlasticDataPtr->mGradPtr = boost::make_shared<MatrixDouble>();
375 commonPlasticDataPtr->mStrainPtr = boost::make_shared<MatrixDouble>();
376 commonPlasticDataPtr->mStressPtr = boost::make_shared<MatrixDouble>();
377
378 CHKERR get_command_line_parameters();
379 CHKERR set_matrial_stiffness();
380
381 if (is_large_strains) {
382 commonHenckyDataPtr = boost::make_shared<HenckyOps::CommonData>();
383 commonHenckyDataPtr->matGradPtr = commonPlasticDataPtr->mGradPtr;
385 commonHenckyDataPtr->matLogCPlastic =
386 commonPlasticDataPtr->getPlasticStrainPtr();
387 commonPlasticDataPtr->mStrainPtr = commonHenckyDataPtr->getMatLogC();
388 commonPlasticDataPtr->mStressPtr =
389 commonHenckyDataPtr->getMatHenckyStress();
390 }
391
393}
constexpr double third
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
double Qinf
Definition: plastic.cpp:108
double rho
Definition: plastic.cpp:101
double visH
Definition: plastic.cpp:104
PetscBool set_timer
Definition: plastic.cpp:95
double scale
Definition: plastic.cpp:97
double zeta
Definition: plastic.cpp:107
double H
Definition: plastic.cpp:103
double cn0
Definition: plastic.cpp:105
double b_iso
Definition: plastic.cpp:109
PetscBool is_large_strains
Definition: plastic.cpp:94
int geom_order
Order if fixed.
Definition: plastic.cpp:112
double sigmaY
Definition: plastic.cpp:102
double cn1
Definition: plastic.cpp:106
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
boost::shared_ptr< PlasticOps::CommonData > commonPlasticDataPtr
Definition: plastic.cpp:156
boost::shared_ptr< HenckyOps::CommonData > commonHenckyDataPtr
Definition: plastic.cpp:157

◆ createCommonData() [2/10]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [3/10]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [4/10]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [5/10]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [6/10]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [7/10]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [8/10]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [9/10]

MoFEMErrorCode Example::createCommonData ( )
private

◆ createCommonData() [10/10]

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}
static Index< 'p', 3 > p
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
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
Definition: integration.cpp:42

◆ 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 = smartCreateDMVector(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,...)
Definition: LogManager.hpp:311
SmartPetscObj< TS > createTSIM(SmartPetscObj< DM > dm=nullptr)
Create TS (time) implicit solver.
const double T

◆ 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]

Examples
plastic.cpp.

Definition at line 459 of file plastic.cpp.

459 {
461 auto pip = mField.getInterface<PipelineManager>();
463 auto bc_mng = mField.getInterface<BcManager>();
464
465 auto add_domain_base_ops = [&](auto &pipeline) {
467
469 "GEOMETRY");
470
471 pipeline.push_back(new OpCalculateScalarFieldValuesDot(
472 "TAU", commonPlasticDataPtr->getPlasticTauDotPtr()));
473 pipeline.push_back(new OpCalculateScalarFieldValues(
474 "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
476 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
478 "EP", commonPlasticDataPtr->getPlasticStrainDotPtr()));
480 "U", commonPlasticDataPtr->mGradPtr));
481
483 };
484
485 auto add_domain_stress_ops = [&](auto &pipeline, auto m_D_ptr) {
487
488 if (is_large_strains) {
489
490 if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
491 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
492 "Wrong pointer for grad");
493
494 pipeline.push_back(
496 pipeline.push_back(
498 pipeline.push_back(
500 pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
501 "U", commonHenckyDataPtr, m_D_ptr));
502 pipeline.push_back(
504
505 } else {
506 pipeline.push_back(
508 commonPlasticDataPtr->mStrainPtr));
509 pipeline.push_back(
510 new OpPlasticStress("U", commonPlasticDataPtr, m_D_ptr, 1));
511 }
512
513 if (m_D_ptr != commonPlasticDataPtr->mDPtr_Axiator) {
514 pipeline.push_back(
516 pipeline.push_back(new OpCalculatePlasticity(
517 "U", mField, commonPlasticDataPtr, m_D_ptr));
518 }
519
521 };
522
523 auto add_domain_ops_lhs_mechanical = [&](auto &pipeline, auto m_D_ptr) {
525 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
526
527 if (is_large_strains) {
528 pipeline.push_back(
530 pipeline.push_back(
531 new OpKPiola("U", "U", commonHenckyDataPtr->getMatTangent()));
533 "U", "EP", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
534 } else {
535 pipeline.push_back(new OpKCauchy("U", "U", m_D_ptr));
536 pipeline.push_back(new OpCalculatePlasticInternalForceLhs_dEP(
537 "U", "EP", commonPlasticDataPtr, m_D_ptr));
538 }
539
540 pipeline.push_back(new OpUnSetBc("U"));
542 };
543
544 auto add_domain_ops_rhs_mechanical = [&](auto &pipeline) {
546 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
547
549 pipeline, mField, "U", {boost::make_shared<PlasticityTimeScale>()},
550 "BODY_FORCE", Sev::inform);
551
552 // Calculate internal forces
553 if (is_large_strains) {
554 pipeline.push_back(new OpInternalForcePiola(
555 "U", commonHenckyDataPtr->getMatFirstPiolaStress()));
556 } else {
557 pipeline.push_back(
558 new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
559 }
560
561 pipeline.push_back(new OpUnSetBc("U"));
563 };
564
565 auto add_domain_ops_lhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
567 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
568
569 if (is_large_strains) {
570 pipeline.push_back(new OpCalculateConstraintsLhs_LogStrain_dU(
571 "TAU", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
572 pipeline.push_back(new OpCalculatePlasticFlowLhs_LogStrain_dU(
573 "EP", "U", commonPlasticDataPtr, commonHenckyDataPtr, m_D_ptr));
574 } else {
575 pipeline.push_back(new OpCalculateConstraintsLhs_dU(
576 "TAU", "U", commonPlasticDataPtr, m_D_ptr));
577 pipeline.push_back(new OpCalculatePlasticFlowLhs_dU(
578 "EP", "U", commonPlasticDataPtr, m_D_ptr));
579 }
580
581 pipeline.push_back(new OpCalculatePlasticFlowLhs_dEP(
582 "EP", "EP", commonPlasticDataPtr, m_D_ptr));
583 pipeline.push_back(new OpCalculatePlasticFlowLhs_dTAU(
584 "EP", "TAU", commonPlasticDataPtr, m_D_ptr));
585 pipeline.push_back(new OpCalculateConstraintsLhs_dEP(
586 "TAU", "EP", commonPlasticDataPtr, m_D_ptr));
587 pipeline.push_back(
589
590 pipeline.push_back(new OpUnSetBc("U"));
592 };
593
594 auto add_domain_ops_rhs_constrain = [&](auto &pipeline, auto m_D_ptr) {
596 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
597
598 pipeline.push_back(
600 pipeline.push_back(
602
603 pipeline.push_back(new OpUnSetBc("U"));
605 };
606
607 auto add_boundary_ops_lhs_mechanical = [&](auto &pipeline) {
611 mField, pipeline, simple->getProblemName(), "U");
613 };
614
615 auto add_boundary_ops_rhs_mechanical = [&](auto &pipeline) {
617
619 pipeline, {NOSPACE}, "GEOMETRY");
620
621 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
623 pipeline, mField, "U", {boost::make_shared<PlasticityTimeScale>()},
624 "FORCE", Sev::inform);
625
626 pipeline.push_back(new OpUnSetBc("U"));
627
628 auto u_mat_ptr = boost::make_shared<MatrixDouble>();
629 pipeline.push_back(
630 new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_mat_ptr));
631
634 mField, pipeline, simple->getProblemName(), "U", u_mat_ptr,
635 {boost::make_shared<TimeScale>()}); // note displacements have no
636 // scaling
637
639 };
640
641 CHKERR add_domain_base_ops(pip->getOpDomainLhsPipeline());
642 CHKERR add_domain_stress_ops(pip->getOpDomainLhsPipeline(),
643 commonPlasticDataPtr->mDPtr);
644 CHKERR add_domain_ops_lhs_mechanical(pip->getOpDomainLhsPipeline(),
645 commonPlasticDataPtr->mDPtr);
646 CHKERR add_domain_ops_lhs_constrain(pip->getOpDomainLhsPipeline(),
647 commonPlasticDataPtr->mDPtr);
648 CHKERR
649 add_boundary_ops_lhs_mechanical(pip->getOpBoundaryLhsPipeline());
650
651 CHKERR add_domain_base_ops(pip->getOpDomainRhsPipeline());
652 CHKERR add_domain_stress_ops(pip->getOpDomainRhsPipeline(),
653 commonPlasticDataPtr->mDPtr);
654 CHKERR add_domain_ops_rhs_mechanical(pip->getOpDomainRhsPipeline());
655 CHKERR add_domain_ops_rhs_constrain(pip->getOpDomainRhsPipeline(),
656 commonPlasticDataPtr->mDPtr);
657
658 // Boundary
659 CHKERR
660 add_boundary_ops_rhs_mechanical(pip->getOpBoundaryRhsPipeline());
661
662 auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
663
664 auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
665
666 CHKERR pip->setDomainRhsIntegrationRule(vol_rule);
667 CHKERR pip->setDomainLhsIntegrationRule(vol_rule);
668
669 CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
670 CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
671
672 auto create_reaction_pipeline = [&](auto &pipeline) {
674
675 if (reactionMarker) {
676
678 "GEOMETRY");
679
680 pipeline.push_back(
682 "U", commonPlasticDataPtr->mGradPtr));
684 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
685
686 if (is_large_strains) {
687
688 if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
689 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
690 "Wrong pointer for grad");
691
692 pipeline.push_back(
694 pipeline.push_back(
696 pipeline.push_back(
698 pipeline.push_back(new OpCalculateHenckyPlasticStress<SPACE_DIM>(
700 pipeline.push_back(
702
703 } else {
704 pipeline.push_back(new OpSymmetrizeTensor<SPACE_DIM>(
705 "U", commonPlasticDataPtr->mGradPtr,
706 commonPlasticDataPtr->mStrainPtr));
707 pipeline.push_back(new OpPlasticStress(
709 }
710
711 pipeline.push_back(new OpSetBc("U", false, reactionMarker));
712 // Calculate internal force
713 if (is_large_strains) {
714 pipeline.push_back(new OpInternalForcePiola(
715 "U", commonHenckyDataPtr->getMatFirstPiolaStress()));
716 } else {
717 pipeline.push_back(
718 new OpInternalForceCauchy("U", commonPlasticDataPtr->mStressPtr));
719 }
720 pipeline.push_back(new OpUnSetBc("U"));
721 }
722
724 };
725
726 reactionFe = boost::make_shared<DomainEle>(mField);
727 reactionFe->getRuleHook = vol_rule;
728
729 CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
730
732}
@ NOSPACE
Definition: definitions.h:83
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForcePiola
Definition: plastic.cpp:68
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradSymTensorGrad< 1, SPACE_DIM, SPACE_DIM, 0 > OpKCauchy
[Only used with Hooke equation (linear material model)]
Definition: plastic.cpp:52
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< G >::OpGradTimesSymTensor< 1, SPACE_DIM, SPACE_DIM > OpInternalForceCauchy
Definition: plastic.cpp:54
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< GAUSS >::OpGradTensorGrad< 1, SPACE_DIM, SPACE_DIM, 1 > OpKPiola
[Only used for dynamics]
Definition: plastic.cpp:66
boost::shared_ptr< DomainEle > reactionFe
Definition: plastic.cpp:158
Add operators pushing bases from local to physical configuration.
Essential boundary conditions.
Definition: Essential.hpp:101
Get value at integration points for scalar field.
Calculate symmetric tensor field rates ant integratio pts.
Calculate symmetric tensor field values at integration pts.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get values at integration pts for tensor filed rank 1, i.e. vector field.

◆ OPs() [2/2]

MoFEMErrorCode Example::OPs ( )
private

◆ outputResults() [1/10]

MoFEMErrorCode Example::outputResults ( )
private

[Solve]

[Postprocess results]

[Postprocessing results]

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

Definition at line 255 of file helmholtz.cpp.

255 {
258 pipeline_mng->getDomainLhsFE().reset();
259 pipeline_mng->getDomainRhsFE().reset();
260 pipeline_mng->getBoundaryLhsFE().reset();
261 pipeline_mng->getBoundaryRhsFE().reset();
262
264
265 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
266
267 auto p_real_ptr = boost::make_shared<VectorDouble>();
268 auto p_imag_ptr = boost::make_shared<VectorDouble>();
269
270 post_proc_fe->getOpPtrVector().push_back(
271 new OpCalculateScalarFieldValues("P_REAL", p_real_ptr));
272 post_proc_fe->getOpPtrVector().push_back(
273 new OpCalculateScalarFieldValues("P_IMAG", p_imag_ptr));
274
276
277 post_proc_fe->getOpPtrVector().push_back(
278
279 new OpPPMap(
280
281 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
282
283 {{"P_REAL", p_real_ptr}, {"P_IMAG", p_imag_ptr}},
284
285 {}, {}, {}
286
287 )
288
289 );
290
291 pipeline_mng->getDomainRhsFE() = post_proc_fe;
292 CHKERR pipeline_mng->loopFiniteElements();
293 CHKERR post_proc_fe->writeFile("out_helmholtz.h5m");
295}
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Post post-proc data at points from hash maps.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()

◆ 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}
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat

◆ 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

◆ 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/10]

MoFEMErrorCode Example::readMesh ( )
private

[run problem]

[Run problem]

[Run programme]

[Read mesh]

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

Definition at line 73 of file helmholtz.cpp.

73 {
75
79
81}
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194

◆ readMesh() [2/10]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [3/10]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [4/10]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [5/10]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [6/10]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [7/10]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [8/10]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [9/10]

MoFEMErrorCode Example::readMesh ( )
private

◆ readMesh() [10/10]

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)
double w(const double g, const double t)
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/13]

MoFEMErrorCode Example::runProblem ( )

[Run problem]

[Create common data]

[Run programme]

[Operator]

[run problem]

[Run all]

[Run problem]

Examples
eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, phase.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 176 of file plastic.cpp.

176 {
180 CHKERR bC();
181 CHKERR OPs();
182 CHKERR tsSolve();
184}
MoFEMErrorCode tsSolve()
Definition: plastic.cpp:752
MoFEMErrorCode createCommonData()
[Set up problem]
Definition: plastic.cpp:261
MoFEMErrorCode OPs()
[Boundary condition]
Definition: plastic.cpp:459
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:188
MoFEMErrorCode bC()
[Create common data]
Definition: plastic.cpp:397

◆ runProblem() [2/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [3/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [4/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [5/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [6/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [7/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [8/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [9/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [10/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [11/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [12/13]

MoFEMErrorCode Example::runProblem ( )

◆ runProblem() [13/13]

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}
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:115
Basic algebra on fields.
Definition: FieldBlas.hpp:21
MoFEMErrorCode setVertexDofs(VertexCoordsFunction lambda, const std::string field_name, Range *verts=nullptr)
Set DOFs on vertices using user function.
Definition: FieldBlas.cpp:320

◆ setIntegrationRules() [1/3]

MoFEMErrorCode Example::setIntegrationRules ( )
private

[Set up problem]

[Set integration rule]

Examples
heat_method.cpp, and plot_base.cpp.

Definition at line 203 of file plot_base.cpp.

203 {
206}

◆ 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 {
140 CHKERR simple->getOptions();
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 .
Definition: definitions.h:60

◆ setupProblem() [1/12]

MoFEMErrorCode Example::setupProblem ( )
private

[Run problem]

[Read mesh]

[Set up problem]

Examples
eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, phase.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 188 of file plastic.cpp.

188 {
191
192 Range domain_ents;
193 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
194 true);
195 auto get_ents_by_dim = [&](const auto dim) {
196 if (dim == SPACE_DIM) {
197 return domain_ents;
198 } else {
199 Range ents;
200 if (dim == 0)
201 CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
202 else
203 CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
204 return ents;
205 }
206 };
207
208 auto get_base = [&]() {
209 auto domain_ents = get_ents_by_dim(SPACE_DIM);
210 if (domain_ents.empty())
211 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
212 const auto type = type_from_handle(domain_ents[0]);
213 switch (type) {
214 case MBQUAD:
216 case MBHEX:
218 case MBTRI:
220 case MBTET:
222 default:
223 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
224 }
225 return NOBASE;
226 };
227
228 const auto base = get_base();
229 MOFEM_LOG("WORLD", Sev::inform) << "Base " << ApproximationBaseNames[base];
230
231 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
232 CHKERR simple->addDomainField("EP", L2, base, size_symm);
233 CHKERR simple->addDomainField("TAU", L2, base, 1);
234 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
235
236 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
237
238 auto ents = get_ents_by_dim(0);
239 ents.merge(get_ents_by_dim(1));
240 // ents.merge(get_ents_by_dim(2));
241 CHKERR simple->setFieldOrder("U", order, &ents);
242 CHKERR simple->setFieldOrder("EP", order - 1);
243 CHKERR simple->setFieldOrder("TAU", order - 2);
244
245 CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
246
247 CHKERR simple->setUp();
248 CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
249
250 auto project_ho_geometry = [&]() {
251 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
252 return mField.loop_dofs("GEOMETRY", ent_method);
253 };
254 CHKERR project_ho_geometry();
255
257}
@ NOBASE
Definition: definitions.h:59
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
static const char *const ApproximationBaseNames[]
Definition: definitions.h:72
const int dim
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
Definition: Templates.hpp:1634
FieldApproximationBase base
Definition: plot_base.cpp:68
Projection of edge entities with one mid-node on hierarchical basis.

◆ setupProblem() [2/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [3/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [4/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [5/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [6/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [7/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [8/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [9/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [10/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [11/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ setupProblem() [12/12]

MoFEMErrorCode Example::setupProblem ( )
private

◆ solveSystem() [1/10]

MoFEMErrorCode Example::solveSystem ( )
private

[Push operators to pipeline]

[Solve]

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

Definition at line 235 of file helmholtz.cpp.

235 {
238 auto solver = pipeline_mng->createKSP();
239 CHKERR KSPSetFromOptions(solver);
240 CHKERR KSPSetUp(solver);
241
242 auto dm = simpleInterface->getDM();
243 auto D = smartCreateDMVector(dm);
244 auto F = smartVectorDuplicate(D);
245
246 CHKERR KSPSolve(solver, F, D);
247 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
248 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
249 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
251}
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.

◆ 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

◆ tsSolve()

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

Definition at line 752 of file plastic.cpp.

752 {
754
757 ISManager *is_manager = mField.getInterface<ISManager>();
758
759 auto snes_ctx_ptr = smartGetDMSnesCtx(simple->getDM());
760
761 auto set_section_monitor = [&](auto solver) {
763 SNES snes;
764 CHKERR TSGetSNES(solver, &snes);
765 CHKERR SNESMonitorSet(snes,
766 (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
768 (void *)(snes_ctx_ptr.get()), nullptr);
770 };
771
772 auto create_post_process_element = [&]() {
773 auto pp_fe = boost::make_shared<PostProcEle>(mField);
775
777 pp_fe->getOpPtrVector(), {H1}, "GEOMETRY");
778
779 auto x_ptr = boost::make_shared<MatrixDouble>();
780 pp_fe->getOpPtrVector().push_back(
781 new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
782 auto u_ptr = boost::make_shared<MatrixDouble>();
783 pp_fe->getOpPtrVector().push_back(
785
786 pp_fe->getOpPtrVector().push_back(
788 "U", commonPlasticDataPtr->mGradPtr));
789 pp_fe->getOpPtrVector().push_back(new OpCalculateScalarFieldValues(
790 "TAU", commonPlasticDataPtr->getPlasticTauPtr()));
791 pp_fe->getOpPtrVector().push_back(
793 "EP", commonPlasticDataPtr->getPlasticStrainPtr()));
794
795 if (is_large_strains) {
796
797 if (commonPlasticDataPtr->mGradPtr != commonHenckyDataPtr->matGradPtr)
798 CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
799
800 pp_fe->getOpPtrVector().push_back(
802 pp_fe->getOpPtrVector().push_back(
804 pp_fe->getOpPtrVector().push_back(
806 pp_fe->getOpPtrVector().push_back(
809 pp_fe->getOpPtrVector().push_back(
811
812 pp_fe->getOpPtrVector().push_back(
813
814 new OpPPMap(
815
816 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
817
818 {},
819
820 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
821
822 {{"GRAD", commonPlasticDataPtr->mGradPtr},
823 {"FIRST_PIOLA", commonHenckyDataPtr->getMatFirstPiolaStress()}},
824
825 {}
826
827 )
828
829 );
830
831 } else {
832 pp_fe->getOpPtrVector().push_back(
834 commonPlasticDataPtr->mStrainPtr));
835 pp_fe->getOpPtrVector().push_back(new OpPlasticStress(
837
838 pp_fe->getOpPtrVector().push_back(
839
840 new OpPPMap(
841
842 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
843
844 {},
845
846 {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
847
848 {},
849
850 {{"STRAIN", commonPlasticDataPtr->mStrainPtr},
851 {"STRESS", commonPlasticDataPtr->mStressPtr}}
852
853 )
854
855 );
856 }
857
858 pp_fe->getOpPtrVector().push_back(
860
861 pp_fe->getOpPtrVector().push_back(
862
863 new OpPPMap(
864
865 pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
866
867 {{"PLASTIC_SURFACE", commonPlasticDataPtr->getPlasticSurfacePtr()},
868 {"PLASTIC_MULTIPLIER", commonPlasticDataPtr->getPlasticTauPtr()}},
869
870 {},
871
872 {},
873
874 {{"PLASTIC_STRAIN", commonPlasticDataPtr->getPlasticStrainPtr()},
875 {"PLASTIC_FLOW", commonPlasticDataPtr->getPlasticFlowPtr()}}
876
877 )
878
879 );
880
881 return pp_fe;
882 };
883
884 auto scatter_create = [&](auto D, auto coeff) {
886 CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
887 ROW, "U", coeff, coeff, is);
888 int loc_size;
889 CHKERR ISGetLocalSize(is, &loc_size);
890 Vec v;
891 CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
892 VecScatter scatter;
893 CHKERR VecScatterCreate(D, is, v, PETSC_NULL, &scatter);
894 return std::make_tuple(SmartPetscObj<Vec>(v),
896 };
897
898 auto set_time_monitor = [&](auto dm, auto solver) {
900 boost::shared_ptr<Monitor> monitor_ptr(
901 new Monitor(dm, create_post_process_element(), reactionFe, uXScatter,
903 boost::shared_ptr<ForcesAndSourcesCore> null;
904 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
905 monitor_ptr, null, null);
907 };
908
909 auto set_fieldsplit_preconditioner = [&](auto solver,
910 boost::shared_ptr<SetUpSchur>
911 &schur_ptr) {
913
914 SNES snes;
915 CHKERR TSGetSNES(solver, &snes);
916 KSP ksp;
917 CHKERR SNESGetKSP(snes, &ksp);
918 PC pc;
919 CHKERR KSPGetPC(ksp, &pc);
920 PetscBool is_pcfs = PETSC_FALSE;
921 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
922
923 // Setup fieldsplit (block) solver - optional: yes/no
924 if (is_pcfs == PETSC_TRUE) {
925
926 auto bc_mng = mField.getInterface<BcManager>();
927 auto name_prb = simple->getProblemName();
928
929 auto create_sub_bc_dm = [&](SmartPetscObj<DM> base_dm,
930 SmartPetscObj<DM> &dm_sub,
931 SmartPetscObj<IS> &is_sub,
932 SmartPetscObj<AO> &ao_sub) {
934
935 dm_sub = createSmartDM(mField.get_comm(), "DMMOFEM");
936 CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SUB_BC");
937 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
938 CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
939 CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
940 for (auto f : {"U", "EP", "TAU"}) {
943 }
944 CHKERR DMSetUp(dm_sub);
945
946 CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_X", "U", 0, 0);
947 CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_Y", "U", 1, 1);
948 CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_Z", "U", 2, 2);
949 CHKERR bc_mng->removeBlockDOFsOnEntities("SUB_BC", "FIX_ALL", "U", 0,
950 2);
951
952 auto *prb_ptr = getProblemPtr(dm_sub);
953 if (auto sub_data = prb_ptr->getSubData()) {
954 is_sub = sub_data->getSmartRowIs();
955 ao_sub = sub_data->getSmartRowMap();
956 int is_sub_size;
957 CHKERR ISGetSize(is_sub, &is_sub_size);
958 MOFEM_LOG("EXAMPLE", Sev::inform)
959 << "Field split second block size " << is_sub_size;
960
961 } else {
962 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "No sub data");
963 }
964
966 };
967
968 auto create_sub_u_dm = [&](SmartPetscObj<DM> base_dm,
969 SmartPetscObj<DM> &dm_sub,
970 SmartPetscObj<IS> &is_sub,
971 SmartPetscObj<AO> &ao_sub) {
973 dm_sub = createSmartDM(mField.get_comm(), "DMMOFEM");
974 CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SUB_U");
975 CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
976 CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
977 CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
978 for (auto f : {"U"}) {
981 }
982 CHKERR DMSetUp(dm_sub);
983
984 auto *prb_ptr = getProblemPtr(dm_sub);
985 if (auto sub_data = prb_ptr->getSubData()) {
986 is_sub = sub_data->getSmartRowIs();
987 ao_sub = sub_data->getSmartRowMap();
988 // ISView(is_sub, PETSC_VIEWER_STDOUT_WORLD);
989 int is_sub_size;
990 CHKERR ISGetSize(is_sub, &is_sub_size);
991 MOFEM_LOG("EXAMPLE", Sev::inform)
992 << "Field split second block size " << is_sub_size;
993
994 } else {
995 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "No sub data");
996 }
997
998 CHKERR DMSetUp(dm_sub);
1000 };
1001
1002 auto create_all_bc_is = [&](SmartPetscObj<IS> &is_all_bc) {
1004 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_X", "U", 0, 0);
1005 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Y", "U", 1, 1, is_all_bc);
1006 is_all_bc = bc_mng->getBlockIS(name_prb, "FIX_Z", "U", 2, 2, is_all_bc);
1007 is_all_bc =
1008 bc_mng->getBlockIS(name_prb, "FIX_ALL", "U", 0, 2, is_all_bc);
1009 int is_all_bc_size;
1010 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
1011 MOFEM_LOG("EXAMPLE", Sev::inform)
1012 << "Field split first block size " << is_all_bc_size;
1014 };
1015
1016 SmartPetscObj<IS> is_all_bc;
1017 SmartPetscObj<DM> dm_bc_sub;
1018 SmartPetscObj<IS> is_bc_sub;
1019 SmartPetscObj<AO> ao_bc_sub;
1020
1021 CHKERR create_all_bc_is(is_all_bc);
1022 CHKERR create_sub_bc_dm(simple->getDM(), dm_bc_sub, is_bc_sub, ao_bc_sub);
1023
1024 SmartPetscObj<IS> is_prb;
1025 CHKERR mField.getInterface<ISManager>()->isCreateProblem(
1026 simple->getProblemName(), ROW, is_prb);
1027
1028 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
1029 is_all_bc); // boundary block
1030 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL, is_bc_sub);
1031
1032 if constexpr (A == AssemblyType::SCHUR) {
1033
1034 SmartPetscObj<IS> is_epp;
1035 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1036 "SUB_BC", ROW, "EP", 0, MAX_DOFS_ON_ENTITY, is_epp);
1037 SmartPetscObj<IS> is_tau;
1038 CHKERR mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
1039 "SUB_BC", ROW, "TAU", 0, MAX_DOFS_ON_ENTITY, is_tau);
1040 IS is_union_raw;
1041 CHKERR ISExpand(is_epp, is_tau, &is_union_raw);
1042 SmartPetscObj<IS> is_union(is_union_raw);
1043
1044 SmartPetscObj<DM> dm_u_sub;
1045 SmartPetscObj<IS> is_u_sub;
1046 SmartPetscObj<AO> ao_u_sub;
1047 CHKERR create_sub_u_dm(dm_bc_sub, dm_u_sub, is_u_sub, ao_u_sub);
1048
1049 // Indices has to be map fro very to level, while assembling Schur
1050 // complement.
1051 auto is_up = is_u_sub;
1052 CHKERR AOPetscToApplicationIS(ao_bc_sub, is_up);
1053 auto ao_up = createAOMappingIS(is_u_sub, PETSC_NULL);
1054 schur_ptr =
1055 SetUpSchur::createSetUpSchur(mField, dm_u_sub, is_union, ao_up);
1056 PetscInt n;
1057 KSP *ksps;
1058 CHKERR PCFieldSplitGetSubKSP(pc, &n, &ksps);
1059 CHKERR schur_ptr->setUp(ksps[1]);
1060 }
1061 }
1062
1064 };
1065
1066 auto dm = simple->getDM();
1067 auto D = smartCreateDMVector(dm);
1068 uXScatter = scatter_create(D, 0);
1069 uYScatter = scatter_create(D, 1);
1070 if constexpr (SPACE_DIM == 3)
1071 uZScatter = scatter_create(D, 2);
1072
1073 auto solver = pip->createTSIM();
1074
1075 auto active_pre_lhs = [&]() {
1077 std::fill(commonPlasticDataPtr->activityData.begin(),
1078 commonPlasticDataPtr->activityData.end(), 0);
1080 };
1081
1082 auto active_post_lhs = [&]() {
1084 auto get_iter = [&]() {
1085 SNES snes;
1086 CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
1087 int iter;
1088 CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
1089 "Can not get iter");
1090 return iter;
1091 };
1092
1093 auto iter = get_iter();
1094 if (iter >= 0) {
1095
1096 std::array<int, 5> activity_data;
1097 std::fill(activity_data.begin(), activity_data.end(), 0);
1098 MPI_Allreduce(commonPlasticDataPtr->activityData.data(),
1099 activity_data.data(), activity_data.size(), MPI_INT,
1100 MPI_SUM, mField.get_comm());
1101
1102 int &active_points = activity_data[0];
1103 int &avtive_full_elems = activity_data[1];
1104 int &avtive_elems = activity_data[2];
1105 int &nb_points = activity_data[3];
1106 int &nb_elements = activity_data[4];
1107
1108 if (nb_points) {
1109
1110 double proc_nb_points =
1111 100 * static_cast<double>(active_points) / nb_points;
1112 double proc_nb_active =
1113 100 * static_cast<double>(avtive_elems) / nb_elements;
1114 double proc_nb_full_active = 100;
1115 if (avtive_elems)
1116 proc_nb_full_active =
1117 100 * static_cast<double>(avtive_full_elems) / avtive_elems;
1118
1120 "EXAMPLE", Sev::inform,
1121 "Iter %d nb pts %d nb avtive pts %d (%3.3f\%) nb active elements %d "
1122 "(%3.3f\%) nb full active elems %d (%3.3f\%)",
1123 iter, nb_points, active_points, proc_nb_points, avtive_elems,
1124 proc_nb_active, avtive_full_elems, proc_nb_full_active, iter);
1125 }
1126 }
1127
1129 };
1130
1131 CHKERR TSSetSolution(solver, D);
1132 CHKERR set_section_monitor(solver);
1133 CHKERR set_time_monitor(dm, solver);
1134 CHKERR TSSetSolution(solver, D);
1135 CHKERR TSSetFromOptions(solver);
1136
1137 boost::shared_ptr<SetUpSchur> schur_ptr;
1138 CHKERR set_fieldsplit_preconditioner(solver, schur_ptr);
1139
1140 mField.getInterface<PipelineManager>()->getDomainLhsFE()->preProcessHook =
1141 [&]() {
1143 if (schur_ptr)
1144 CHKERR schur_ptr->preProc();
1145 CHKERR active_pre_lhs();
1147 };
1148 mField.getInterface<PipelineManager>()->getDomainLhsFE()->postProcessHook =
1149 [&]() {
1151 // if (schur_ptr)
1152 // CHKERR schur_ptr->postProc();
1153 CHKERR active_post_lhs();
1155 };
1156
1157 mField.getInterface<PipelineManager>()->getBoundaryLhsFE()->postProcessHook =
1158 [&]() {
1160 if (schur_ptr)
1161 CHKERR schur_ptr->postProc();
1163 };
1164
1165 MOFEM_LOG_CHANNEL("TIMER");
1166 MOFEM_LOG_TAG("TIMER", "timer");
1167 if(set_timer)
1168 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
1169 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
1170 CHKERR TSSetUp(solver);
1171 MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
1172 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
1173 CHKERR TSSolve(solver, NULL);
1174 MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
1175
1176 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
1177 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
1178 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
1179
1181}
#define MAX_DOFS_ON_ENTITY
Maximal number of DOFs on entity.
Definition: definitions.h:236
FTensor::Index< 'n', SPACE_DIM > n
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:221
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:485
PetscErrorCode DMMoFEMSetSquareProblem(DM dm, PetscBool square_problem)
set squared problem
Definition: DMMoFEM.cpp:444
PetscErrorCode DMMoFEMAddSubFieldRow(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:244
PetscErrorCode DMMoFEMAddSubFieldCol(DM dm, const char field_name[])
Definition: DMMoFEM.cpp:277
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:339
const FTensor::Tensor2< T, Dim, Dim > Vec
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:1044
MoFEMErrorCode MoFEMSNESMonitorFields(SNES snes, PetscInt its, PetscReal fgnorm, SnesCtx *snes_ctx)
Sens monitor printing residual field by field.
Definition: SnesCtx.cpp:224
auto createSmartDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
auto createAOMappingIS(IS isapp, IS ispetsc)
Creates an application mapping using two index sets.
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition: DMMoFEM.hpp:971
auto smartGetDMSnesCtx(DM dm)
Get SNES context data structure used by DM.
Definition: DMMoFEM.hpp:1017
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uZScatter
Definition: plastic.cpp:162
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uXScatter
Definition: plastic.cpp:160
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > > uYScatter
Definition: plastic.cpp:161
Section manager is used to create indexes and sections.
Definition: ISManager.hpp:23
intrusive_ptr for managing petsc objects
static boost::shared_ptr< SetUpSchur > createSetUpSchur(MoFEM::Interface &m_field, SmartPetscObj< DM >, SmartPetscObj< IS >, SmartPetscObj< AO >)
Definition: plastic.cpp:1351

Member Data Documentation

◆ 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

FieldApproximationBase Example::base
private
Examples
plastic.cpp, and plot_base.cpp.

Definition at line 68 of file plot_base.cpp.

◆ boundaryMarker

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

Definition at line 164 of file plastic.cpp.

◆ commonDataPtr

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

Definition at line 42 of file integration.cpp.

◆ commonHenckyDataPtr

boost::shared_ptr< HenckyOps::CommonData > Example::commonHenckyDataPtr
private
Examples
nonlinear_elastic.cpp, and plastic.cpp.

Definition at line 157 of file plastic.cpp.

◆ commonPlasticDataPtr

boost::shared_ptr<PlasticOps::CommonData> Example::commonPlasticDataPtr
private
Examples
plastic.cpp.

Definition at line 156 of file plastic.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.

◆ 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.

◆ K

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

Definition at line 68 of file eigen_elastic.cpp.

◆ M

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

Definition at line 67 of file eigen_elastic.cpp.

◆ matAccelerationPtr

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

Definition at line 77 of file nonlinear_dynamic_elastic.cpp.

◆ matDPtr

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

Definition at line 65 of file eigen_elastic.cpp.

◆ matGradPtr

boost::shared_ptr< MatrixDouble > Example::matGradPtr
private
Examples
nonlinear_elastic.cpp.

Definition at line 76 of file nonlinear_elastic.cpp.

◆ matInertiaPtr

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

Definition at line 78 of file nonlinear_dynamic_elastic.cpp.

◆ matStrainPtr

boost::shared_ptr< MatrixDouble > Example::matStrainPtr
private
Examples
nonlinear_elastic.cpp.

Definition at line 77 of file nonlinear_elastic.cpp.

◆ matStressPtr

boost::shared_ptr< MatrixDouble > Example::matStressPtr
private
Examples
nonlinear_elastic.cpp.

Definition at line 78 of file nonlinear_elastic.cpp.

◆ matTangentPtr

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

Definition at line 81 of file nonlinear_dynamic_elastic.cpp.

◆ mField

MoFEM::Interface & Example::mField
private

◆ pinchNodes

Range Example::pinchNodes
private
Examples
heat_method.cpp.

Definition at line 54 of file heat_method.cpp.

◆ reactionFe

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

Definition at line 158 of file plastic.cpp.

◆ reactionMarker

boost::shared_ptr<std::vector<unsigned char> > Example::reactionMarker
private
Examples
plastic.cpp.

Definition at line 165 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 45 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 160 of file plastic.cpp.

◆ uYScatter

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

Definition at line 161 of file plastic.cpp.

◆ uZScatter

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

Definition at line 162 of file plastic.cpp.


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