v0.16.3
Loading...
Searching...
No Matches
Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | List of all members
ElasticExample Struct Reference

#include "tutorials/vec-0_elasticity/src/ElasticExample.hpp"

Inheritance diagram for ElasticExample:
[legend]
Collaboration diagram for ElasticExample:
[legend]

Public Member Functions

 ElasticExample (MoFEM::Interface &m_field)
 
virtual ~ElasticExample ()=default
 
MoFEMErrorCode runProblem ()
 [Run problem]
 

Protected Types

using Pipeline = boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator >
 

Protected Member Functions

virtual MoFEMErrorCode readMesh ()
 [Run problem]
 
virtual MoFEMErrorCode setupProblem ()
 [Read mesh]
 
virtual MoFEMErrorCode boundaryCondition ()
 [Set up problem]
 
virtual MoFEMErrorCode assembleSystem ()
 
virtual MoFEMErrorCode solveSystem ()
 [Solve]
 
virtual MoFEMErrorCode outputResults ()
 [Solve]
 
virtual MoFEMErrorCode checkResults ()
 [Postprocess results]
 
virtual MoFEMErrorCode addDomainMaterialOps (Pipeline &lhs_pipeline, Pipeline &rhs_pipeline)
 Add constitutive operators to the domain assembly pipelines.
 
virtual MoFEMErrorCode setUpSolver (SmartPetscObj< KSP > solver)
 [Push operators to pipeline]
 
virtual MoFEMErrorCode kspSetUpAndSolve (SmartPetscObj< KSP > solver)
 

Protected Attributes

MoFEM::Interface & mField
 
boost::shared_ptr< MatrixDouble > vectorFieldPtr = nullptr
 

Detailed Description

Definition at line 18 of file ElasticExample.hpp.

Member Typedef Documentation

◆ Pipeline

using ElasticExample::Pipeline = boost::ptr_deque<ForcesAndSourcesCore::UserDataOperator>
protected

Definition at line 26 of file ElasticExample.hpp.

Constructor & Destructor Documentation

◆ ElasticExample()

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

Definition at line 20 of file ElasticExample.hpp.

20: mField(m_field) {}

◆ ~ElasticExample()

virtual ElasticExample::~ElasticExample ( )
virtualdefault

Member Function Documentation

◆ addDomainMaterialOps()

MoFEMErrorCode ElasticExample::addDomainMaterialOps ( Pipeline &  lhs_pipeline,
Pipeline &  rhs_pipeline 
)
protectedvirtual

Add constitutive operators to the domain assembly pipelines.

[Boundary condition]

[Push operators to pipeline]

Definition at line 147 of file ElasticExample.hpp.

148 {
150 // Add LHS operator for elasticity (stiffness matrix)
151 CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
152 mField, lhs_pipeline, "U", "MAT_ELASTIC", Sev::verbose);
153 // Add RHS operator for internal forces
154 CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
155 mField, rhs_pipeline, "U", "MAT_ELASTIC", Sev::verbose);
157}
#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.

◆ assembleSystem()

MoFEMErrorCode ElasticExample::assembleSystem ( )
protectedvirtual

[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
mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 159 of file ElasticExample.hpp.

159 {
161 auto pip = mField.getInterface<PipelineManager>();
162
163 //! [Integration rule]
164 auto integration_rule = [](int, int, int approx_order) {
165 return 2 * approx_order + 1;
166 };
167
168 auto integration_rule_bc = [](int, int, int approx_order) {
169 return 2 * approx_order + 1;
170 };
171
172 CHKERR pip->setDomainRhsIntegrationRule(integration_rule);
173 CHKERR pip->setDomainLhsIntegrationRule(integration_rule);
174 CHKERR pip->setBoundaryRhsIntegrationRule(integration_rule_bc);
175 CHKERR pip->setBoundaryLhsIntegrationRule(integration_rule_bc);
176 //! [Integration rule]
177
179 pip->getOpDomainLhsPipeline(), {H1}, "GEOMETRY");
181 pip->getOpDomainRhsPipeline(), {H1}, "GEOMETRY");
183 pip->getOpBoundaryRhsPipeline(), {NOSPACE}, "GEOMETRY");
185 pip->getOpBoundaryLhsPipeline(), {NOSPACE}, "GEOMETRY");
186
187 //! [Push domain stiffness matrix to pipeline]
188 CHKERR addDomainMaterialOps(pip->getOpDomainLhsPipeline(),
189 pip->getOpDomainRhsPipeline());
190 //! [Push domain stiffness matrix to pipeline]
191
192 //! [Push Body forces]
193 CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
194 pip->getOpDomainRhsPipeline(), mField, "U", Sev::inform);
195 //! [Push Body forces]
196
197 //! [Push natural boundary conditions]
198 // Add force boundary condition
199 CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
200 pip->getOpBoundaryRhsPipeline(), mField, "U", -1, Sev::inform);
201 // Add case for mix type of BCs
202 CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
203 pip->getOpBoundaryLhsPipeline(), mField, "U", Sev::verbose);
204 //! [Push natural boundary conditions]
206}
auto integration_rule
static constexpr int approx_order
virtual MoFEMErrorCode addDomainMaterialOps(Pipeline &lhs_pipeline, Pipeline &rhs_pipeline)
Add constitutive operators to the domain assembly pipelines.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.

◆ boundaryCondition()

MoFEMErrorCode ElasticExample::boundaryCondition ( )
protectedvirtual

[Set up problem]

[Boundary condition]

Examples
mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 125 of file ElasticExample.hpp.

125 {
127 auto simple = mField.getInterface<Simple>();
128 auto bc_mng = mField.getInterface<BcManager>();
129
130 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
131 "U", 0, 0);
132 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
133 "U", 1, 1);
134 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
135 "U", 2, 2);
136 CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
137 "REMOVE_ALL", "U", 0, 3);
138 CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
139 simple->getProblemName(), "U");
140
142}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition acoustic.cpp:69

◆ checkResults()

MoFEMErrorCode ElasticExample::checkResults ( )
protectedvirtual

[Postprocess results]

[Check]

Reimplemented in ElasticTieMeshTutorial, and ElasticTieMeshExample.

Examples
mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 362 of file ElasticExample.hpp.

362 {
363 MOFEM_LOG_CHANNEL("WORLD");
364 auto simple = mField.getInterface<Simple>();
365 auto pip = mField.getInterface<PipelineManager>();
367
368 auto integration_rule = [](int, int, int p_data) { return 2 * p_data + 1; };
369 CHKERR pip->setEvaluationIntegrationRule(integration_rule);
370
371 auto &domain_evaluation_pipeline = pip->getOpEvaluationPipeline();
372
373 if (domain_evaluation_pipeline.empty()) {
375 domain_evaluation_pipeline, {H1}, "GEOMETRY");
376
377 auto op_loop_boundary = new OpLoopSide<BoundaryEle>(
378 mField, simple->getBoundaryFEName(), SPACE_DIM - 1, Sev::verbose);
379 op_loop_boundary->getSideFEPtr()->getRuleHook = integration_rule;
380 auto &boundary_evaluation_pipeline = op_loop_boundary->getOpPtrVector();
381
383 boundary_evaluation_pipeline, {}, "GEOMETRY");
384
385 // Add RHS operators for internal forces
386 CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
387 mField, domain_evaluation_pipeline, "U", "MAT_ELASTIC", Sev::verbose);
388
389 CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
390 domain_evaluation_pipeline, mField, "U", Sev::verbose);
391
392 CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
393 boundary_evaluation_pipeline, mField, "U", -1, Sev::verbose);
394 domain_evaluation_pipeline.push_back(op_loop_boundary);
395 }
396
397 auto dm = simple->getDM();
398 auto res = createDMVector(dm);
399 CHKERR VecSetDM(res, PETSC_NULLPTR);
400
401 pip->getEvaluationFE()->f = res;
402
403 CHKERR VecZeroEntries(res);
404
405 CHKERR pip->loopFiniteElementsEvaluation();
406 // CHKERR mField.getInterface<FieldBlas>()->fieldScale(-1, "U");
407
408 CHKERR VecGhostUpdateBegin(res, ADD_VALUES, SCATTER_REVERSE);
409 CHKERR VecGhostUpdateEnd(res, ADD_VALUES, SCATTER_REVERSE);
410 CHKERR VecAssemblyBegin(res);
411 CHKERR VecAssemblyEnd(res);
412
413 auto zero_residual_at_constrains = [&]() {
415 auto fe_post_proc_ptr = boost::make_shared<FEMethod>();
416 auto get_post_proc_hook_rhs = [&]() {
418 CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
419 mField, fe_post_proc_ptr, res)();
420 CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
421 mField, fe_post_proc_ptr, 0, res)();
423 };
424 fe_post_proc_ptr->postProcessHook = get_post_proc_hook_rhs;
425 CHKERR DMoFEMPostProcessFiniteElements(dm, fe_post_proc_ptr.get());
427 };
428
429 CHKERR zero_residual_at_constrains();
430
431 double nrm2;
432 CHKERR VecNorm(res, NORM_2, &nrm2);
433 MOFEM_LOG_CHANNEL("WORLD");
434 MOFEM_LOG_C("WORLD", Sev::inform, "residual = %3.4e\n", nrm2);
435
436 int test = 0;
437 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-test", &test, PETSC_NULLPTR);
438 if (test > 0) {
439 auto post_proc_residual = [&](auto dm, auto f_res, auto out_name) {
441 auto post_proc_fe =
442 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
443 using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
444 auto u_vec = boost::make_shared<MatrixDouble>();
445 post_proc_fe->getOpPtrVector().push_back(
446 new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_vec, f_res));
447 post_proc_fe->getOpPtrVector().push_back(
448
449 new OpPPMap(
450
451 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
452
453 {},
454
455 {{"RES", u_vec}},
456
457 {}, {})
458
459 );
460
461 CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
462 post_proc_fe);
463 post_proc_fe->writeFile(out_name);
465 };
466
467 CHKERR post_proc_residual(simple->getDM(), res, "res.h5m");
468
469 constexpr double eps = 1e-8;
470 if (nrm2 > eps)
471 SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
472 "Residual is not zero");
473 }
474 if (test == 2) {
475 if (!vectorFieldPtr || vectorFieldPtr->size1() == 0) {
476 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
477 "atom test %d failed: Field Evaluator did not provide result",
478 test);
479 }
480 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
481 double Ux_ref = 0.46;
482 double Uy_ref = -0.03;
483 constexpr double eps = 1e-8;
484 if (fabs(t_disp(0) - Ux_ref) > eps || fabs(t_disp(1) - Uy_ref) > eps) {
485 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
486 "atom test %d failed: Ux_ref = %3.6e, computed = %3.6e, Uy_ref "
487 "= %3.6e, computed = %3.6e",
488 test, Ux_ref, t_disp(0), Uy_ref, t_disp(1));
489 }
490 }
492}
#define MOFEM_LOG_C(channel, severity, format,...)
static const double eps
constexpr int SPACE_DIM
@ H1
continuous field
Definition definitions.h:85
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
PetscErrorCode DMoFEMPostProcessFiniteElements(DM dm, MoFEM::FEMethod *method)
execute finite element method for each element in dm (problem)
Definition DMMoFEM.cpp:546
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, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1237
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
MoFEM::Interface & mField
boost::shared_ptr< MatrixDouble > vectorFieldPtr

◆ kspSetUpAndSolve()

MoFEMErrorCode ElasticExample::kspSetUpAndSolve ( SmartPetscObj< KSP >  solver)
protectedvirtual

Reimplemented in ElasticSchurExample, ElasticTieMeshExample, and ElasticAdaptiveExample.

Definition at line 212 of file ElasticExample.hpp.

212 {
214
215 MOFEM_LOG_CHANNEL("TIMER");
216 MOFEM_LOG_TAG("TIMER", "timer");
217
218 BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
219 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp";
220 CHKERR KSPSetUp(solver);
221 MOFEM_LOG("TIMER", Sev::inform) << "KSPSetUp <= Done";
222
223 DM dm;
224 CHKERR KSPGetDM(solver, &dm);
225 auto D = createDMVector(dm);
226 auto F = vectorDuplicate(D);
227
228 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve";
229 CHKERR KSPSolve(solver, F, D);
230 MOFEM_LOG("TIMER", Sev::inform) << "KSPSolve <= Done";
231
232 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
233 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
234 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
235
237};
@ F
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
double D

◆ outputResults()

MoFEMErrorCode ElasticExample::outputResults ( )
protectedvirtual

[Solve]

[Postprocess results]

Examples
mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 351 of file ElasticExample.hpp.

351 {
353 auto simple = mField.getInterface<Simple>();
356 mField, simple->getDM(), simple->getDomainFEName(), "out_elastic.h5m");
358}
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
MoFEMErrorCode postProcessElasticResults(MoFEM::Interface &mField, SmartPetscObj< DM > dm, const std::string &domain_fe_name, const std::string &out_file_name, std::vector< std::pair< std::string, SmartPetscObj< Vec > > > extra_vectors={}, const std::vector< std::string > &tags_to_transfer={}, const Sev hooke_ops_sev=Sev::verbose)
ElementsAndOps< SPACE_DIM >::SideEle SideEle
Definition plastic.cpp:61

◆ readMesh()

MoFEMErrorCode ElasticExample::readMesh ( )
protectedvirtual

[Run problem]

[Read mesh]

Reimplemented in ElasticAdaptiveExample, and ElasticTieMeshExample.

Examples
mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 63 of file ElasticExample.hpp.

63 {
65 auto simple = mField.getInterface<Simple>();
66 CHKERR simple->getOptions();
67 CHKERR simple->loadFile();
68 // Add meshsets if config file provided
69 CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
71}

◆ runProblem()

MoFEMErrorCode ElasticExample::runProblem ( )

[Run problem]

Definition at line 49 of file ElasticExample.hpp.

49 {
51 CHKERR readMesh();
52 CHKERR setupProblem();
53 CHKERR boundaryCondition();
54 CHKERR assembleSystem();
55 CHKERR solveSystem();
56 CHKERR outputResults();
57 CHKERR checkResults();
59}

◆ setupProblem()

MoFEMErrorCode ElasticExample::setupProblem ( )
protectedvirtual

[Read mesh]

[Set up problem]

Reimplemented in ElasticTieMeshExample.

Examples
mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 75 of file ElasticExample.hpp.

75 {
77 Simple *simple = mField.getInterface<Simple>();
78
79 enum bases { AINSWORTH, DEMKOWICZ, LASBASETOPT };
80 const char *list_bases[LASBASETOPT] = {"ainsworth", "demkowicz"};
81 PetscInt choice_base_value = AINSWORTH;
82 CHKERR PetscOptionsGetEList(PETSC_NULLPTR, NULL, "-base", list_bases,
83 LASBASETOPT, &choice_base_value, PETSC_NULLPTR);
84
86 switch (choice_base_value) {
87 case AINSWORTH:
89 MOFEM_LOG("WORLD", Sev::inform)
90 << "Set AINSWORTH_LEGENDRE_BASE for displacements";
91 break;
92 case DEMKOWICZ:
94 MOFEM_LOG("WORLD", Sev::inform)
95 << "Set DEMKOWICZ_JACOBI_BASE for displacements";
96 break;
97 default:
98 base = LASTBASE;
99 break;
100 }
101
102 // Add field
103 CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
104 CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
105 int order = 2;
106 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
107
108 CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
109
110 CHKERR simple->setFieldOrder("U", order);
111 CHKERR simple->setFieldOrder("GEOMETRY", 2);
112 CHKERR simple->setUp();
113
114 auto project_ho_geometry = [&]() {
115 Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
116 return mField.loop_dofs("GEOMETRY", ent_method);
117 };
118 CHKERR project_ho_geometry();
119
121}
FieldApproximationBase
approximation base
Definition definitions.h:58
@ LASTBASE
Definition definitions.h:69
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
constexpr int order
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.
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)

◆ setUpSolver()

MoFEMErrorCode ElasticExample::setUpSolver ( SmartPetscObj< KSP >  solver)
protectedvirtual

[Push operators to pipeline]

Reimplemented in ElasticTieMeshExample.

Definition at line 209 of file ElasticExample.hpp.

209{ return 0; }

◆ solveSystem()

MoFEMErrorCode ElasticExample::solveSystem ( )
protectedvirtual

[Solve]

Examples
mofem/tutorials/vec-10_elasticity_schur/schur_elastic.cpp.

Definition at line 241 of file ElasticExample.hpp.

241 {
243
244 auto simple = mField.getInterface<Simple>();
245 auto pip = mField.getInterface<PipelineManager>();
246 auto solver = pip->createKSP();
247 CHKERR setUpSolver(solver);
248 CHKERR KSPSetFromOptions(solver);
249
250 auto set_essential_bc = [this]() {
252 // This is low level pushing finite elements (pipelines) to solver
253 auto simple = mField.getInterface<Simple>();
254 auto dm = simple->getDM();
255 auto ksp_ctx_ptr = getDMKspCtx(dm);
256
257 auto pre_proc_rhs = boost::make_shared<FEMethod>();
258 auto post_proc_rhs = boost::make_shared<FEMethod>();
259 auto post_proc_lhs = boost::make_shared<FEMethod>();
260
261 auto get_pre_proc_hook = [this, pre_proc_rhs]() {
262 return EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_rhs,
263 {});
264 };
265 pre_proc_rhs->preProcessHook = get_pre_proc_hook();
266
267 auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
269
270 CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
271 post_proc_rhs, 1.)();
273 };
274
275 auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
277
278 CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(mField,
279 post_proc_lhs, 1.)();
281 };
282
283 post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
284 post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
285
286 ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
287 ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
288 ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
290 };
291
292 auto evaluate_field_at_the_point = [&]() {
294
295 int coords_dim = 3;
296 std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
297 PetscBool do_eval_field = PETSC_FALSE;
298 CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
299 field_eval_coords.data(), &coords_dim,
301
302 if (do_eval_field) {
303
304 vectorFieldPtr = boost::make_shared<MatrixDouble>();
305 auto field_eval_data =
306 mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
307
308 CHKERR mField.getInterface<FieldEvaluatorInterface>()
309 ->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
310
311 field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
312 auto no_rule = [](int, int, int) { return -1; };
313 auto field_eval_fe_ptr = field_eval_data->feMethodPtr;
314 field_eval_fe_ptr->getRuleHook = no_rule;
315
316 field_eval_fe_ptr->getOpPtrVector().push_back(
317 new OpCalculateVectorFieldValues<SPACE_DIM>("U", vectorFieldPtr));
318
319 CHKERR mField.getInterface<FieldEvaluatorInterface>()
320 ->evalFEAtThePoint<SPACE_DIM>(
321 field_eval_coords.data(), 1e-12, simple->getProblemName(),
322 simple->getDomainFEName(), field_eval_data,
324 QUIET);
325
326 if (vectorFieldPtr->size1()) {
327 auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
328 if constexpr (SPACE_DIM == 2)
329 MOFEM_LOG("FieldEvaluator", Sev::inform)
330 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
331 else
332 MOFEM_LOG("FieldEvaluator", Sev::inform)
333 << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
334 << " U_Z: " << t_disp(2);
335 }
336
338 }
340 };
341
342 CHKERR set_essential_bc();
343 CHKERR kspSetUpAndSolve(solver);
344 CHKERR evaluate_field_at_the_point();
345
347}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
@ QUIET
@ MF_EXIST
#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 PetscOptionsGetRealArray(PetscOptions *, const char pre[], const char name[], PetscReal dval[], PetscInt *nmax, PetscBool *set)
virtual MoFEMErrorCode kspSetUpAndSolve(SmartPetscObj< KSP > solver)
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
PetscBool do_eval_field
Evaluate field.
Definition plastic.cpp:119

Member Data Documentation

◆ mField

MoFEM::Interface& ElasticExample::mField
protected

◆ vectorFieldPtr

boost::shared_ptr<MatrixDouble> ElasticExample::vectorFieldPtr = nullptr
protected

Definition at line 30 of file ElasticExample.hpp.


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