v0.13.1
heat_equation.cpp
Go to the documentation of this file.
1/**
2 * \file heat_equation.cpp
3 * \example heat_equation.cpp
4 *
5 * \brief Solve the time-dependent Heat Equation
6 \f[
7 \begin{aligned}
8\frac{\partial u(\mathbf{x}, t)}{\partial t}-\Delta u(\mathbf{x}, t)
9&=f(\mathbf{x}, t) & & \forall \mathbf{x} \in \Omega, t \in(0, T), \\
10u(\mathbf{x}, 0) &=u_{0}(\mathbf{x}) & & \forall \mathbf{x} \in \Omega, \\
11u(\mathbf{x}, t) &=g(\mathbf{x}, t) & & \forall \mathbf{x} \in \partial \Omega,
12t \in(0, T). \end{aligned}
13 \f]
14 **/
15
16/* This file is part of MoFEM.
17 * MoFEM is free software: you can redistribute it and/or modify it under
18 * the terms of the GNU Lesser General Public License as published by the
19 * Free Software Foundation, either version 3 of the License, or (at your
20 * option) any later version.
21 *
22 * MoFEM is distributed in the hope that it will be useful, but WITHOUT
23 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
24 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
25 * License for more details.
26 *
27 * You should have received a copy of the GNU Lesser General Public
28 * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
29
30#include <stdlib.h>
31#include <cmath>
32#include <BasicFiniteElements.hpp>
33
34using namespace MoFEM;
35
36static char help[] = "...\n\n";
37
38template <int DIM> struct ElementsAndOps {};
39
40//! [Define dimension]
41constexpr int SPACE_DIM = 2; //< Space dimension of problem, mesh
42//! [Define dimension]
43
50
61
68
69// Capacity
70constexpr double c = 1;
71constexpr double k = 1;
72constexpr double init_u = 0.;
73
74/**
75 * @brief Monitor solution
76 *
77 * This functions is called by TS solver at the end of each step. It is used
78 * to output results to the hard drive.
79 */
80struct Monitor : public FEMethod {
81
82 Monitor(SmartPetscObj<DM> dm, boost::shared_ptr<PostProcEle> post_proc)
83 : dM(dm), postProc(post_proc){};
84
85 MoFEMErrorCode preProcess() { return 0; }
86 MoFEMErrorCode operator()() { return 0; }
87
88 static constexpr int saveEveryNthStep = 1;
89
92 if (ts_step % saveEveryNthStep == 0) {
94 CHKERR postProc->writeFile(
95 "out_level_" + boost::lexical_cast<std::string>(ts_step) + ".h5m");
96 }
98 }
99
100private:
102 boost::shared_ptr<PostProcEle> postProc;
103};
104
106public:
108
109 // Declaration of the main function to run analysis
111
112private:
113 // Declaration of other main functions called in runProgram()
122
123 // Main interfaces
125
126 // Object to mark boundary entities for the assembling of domain elements
127 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
128};
129
130HeatEquation::HeatEquation(MoFEM::Interface &m_field) : mField(m_field) {}
131
134
135 auto *simple = mField.getInterface<Simple>();
137 CHKERR simple->getOptions();
138 CHKERR simple->loadFile();
139
141}
142
145
146 auto *simple = mField.getInterface<Simple>();
147 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
148 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
149
150 int order = 3;
151 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
152 CHKERR simple->setFieldOrder("U", order);
153
154 CHKERR simple->setUp();
155
157}
158
161
162 auto integration_rule = [](int o_row, int o_col, int approx_order) {
163 return 2 * approx_order;
164 };
165
166 auto *pipeline_mng = mField.getInterface<PipelineManager>();
167 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
168 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
169 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
170 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
171
173}
174
177
178 // Get surface entities form blockset, set initial values in those
179 // blocksets. To keep it simple, it is assumed that inital values are on
180 // blockset 1
182 Range inner_surface;
183 CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
184 1, BLOCKSET, 2, inner_surface, true);
185 if (!inner_surface.empty()) {
186 Range inner_surface_verts;
187 CHKERR mField.get_moab().get_connectivity(inner_surface,
188 inner_surface_verts, false);
189 CHKERR mField.getInterface<FieldBlas>()->setField(
190 init_u, MBVERTEX, inner_surface_verts, "U");
191 }
192 }
193
195}
196
199
200 auto bc_mng = mField.getInterface<BcManager>();
201 auto *simple = mField.getInterface<Simple>();
202 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ESSENTIAL",
203 "U", 0, 0);
204
205 auto &bc_map = bc_mng->getBcMapByBlockName();
206 boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
207 for (auto b : bc_map) {
208 if (std::regex_match(b.first, std::regex("(.*)ESSENTIAL(.*)"))) {
209 boundaryMarker->resize(b.second->bcMarkers.size(), 0);
210 for (int i = 0; i != b.second->bcMarkers.size(); ++i) {
211 (*boundaryMarker)[i] |= b.second->bcMarkers[i];
212 }
213 }
214 }
215
217}
218
221
222 auto add_domain_base_ops = [&](auto &pipeline) {
223 auto det_ptr = boost::make_shared<VectorDouble>();
224 auto jac_ptr = boost::make_shared<MatrixDouble>();
225 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
226 pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
227 pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
228 pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
229 pipeline.push_back(new OpSetHOWeights(det_ptr));
230 };
231
232 auto add_domain_lhs_ops = [&](auto &pipeline) {
233 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
234 pipeline.push_back(new OpDomainGradGrad(
235 "U", "U", [](double, double, double) -> double { return k; }));
236 auto get_c = [this](const double, const double, const double) {
237 auto pipeline_mng = mField.getInterface<PipelineManager>();
238 auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
239 return c * fe_domain_lhs->ts_a;
240 };
241 pipeline.push_back(new OpDomainMass("U", "U", get_c));
242 pipeline.push_back(new OpUnSetBc("U"));
243 };
244
245 auto add_domain_rhs_ops = [&](auto &pipeline) {
246 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
247 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
248 auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
249 pipeline.push_back(new OpCalculateScalarFieldGradient<SPACE_DIM>(
250 "U", grad_u_at_gauss_pts));
251 pipeline.push_back(
252 new OpCalculateScalarFieldValuesDot("U", dot_u_at_gauss_pts));
253 pipeline.push_back(new OpDomainGradTimesVec(
254 "U", grad_u_at_gauss_pts,
255 [](double, double, double) -> double { return k; }));
256 pipeline.push_back(new OpDomainTimesScalarField(
257 "U", dot_u_at_gauss_pts,
258 [](const double, const double, const double) { return c; }));
259 auto source_term = [&](const double x, const double y, const double z) {
260 auto pipeline_mng = mField.getInterface<PipelineManager>();
261 auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
262 const auto t = fe_domain_lhs->ts_t;
263 return 1e1 * pow(M_E, -M_PI * M_PI * t) * sin(1. * M_PI * x) *
264 sin(2. * M_PI * y);
265 };
266 pipeline.push_back(new OpDomainSource("U", source_term));
267 pipeline.push_back(new OpUnSetBc("U"));
268 };
269
270 auto add_boundary_base_ops = [&](auto &pipeline) {};
271
272 auto add_lhs_base_ops = [&](auto &pipeline) {
273 pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
274 pipeline.push_back(new OpBoundaryMass(
275 "U", "U", [](const double, const double, const double) { return c; }));
276 pipeline.push_back(new OpUnSetBc("U"));
277 };
278 auto add_rhs_base_ops = [&](auto &pipeline) {
279 pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
280 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
281 auto boundary_function = [&](const double x, const double y,
282 const double z) {
283 auto pipeline_mng = mField.getInterface<PipelineManager>();
284 auto &fe_rhs = pipeline_mng->getBoundaryRhsFE();
285 const auto t = fe_rhs->ts_t;
286 return 0;
287 // abs(0.1 * pow(M_E, -M_PI * M_PI * t) * sin(2. * M_PI * x) *
288 // sin(3. * M_PI * y));
289 };
290 pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
291 pipeline.push_back(new OpBoundaryTimeScalarField(
292 "U", u_at_gauss_pts,
293 [](const double, const double, const double) { return c; }));
294 pipeline.push_back(new OpBoundarySource("U", boundary_function));
295 pipeline.push_back(new OpUnSetBc("U"));
296 };
297
298 auto pipeline_mng = mField.getInterface<PipelineManager>();
299 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
300 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
301 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
302 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
303
304 add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
305 add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
306 add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
307 add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
308
310}
311
314
315 auto *simple = mField.getInterface<Simple>();
316 auto *pipeline_mng = mField.getInterface<PipelineManager>();
317
318 auto create_post_process_element = [&]() {
319 auto post_froc_fe = boost::make_shared<PostProcEle>(mField);
320 post_froc_fe->generateReferenceElementMesh();
321 auto det_ptr = boost::make_shared<VectorDouble>();
322 auto jac_ptr = boost::make_shared<MatrixDouble>();
323 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
324 post_froc_fe->getOpPtrVector().push_back(new OpCalculateHOJac<2>(jac_ptr));
325 post_froc_fe->getOpPtrVector().push_back(
326 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
327 post_froc_fe->getOpPtrVector().push_back(
328 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
329 post_froc_fe->addFieldValuesPostProc("U");
330 return post_froc_fe;
331 };
332
333 auto set_time_monitor = [&](auto dm, auto solver) {
335 boost::shared_ptr<Monitor> monitor_ptr(
336 new Monitor(dm, create_post_process_element()));
337 boost::shared_ptr<ForcesAndSourcesCore> null;
338 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
339 monitor_ptr, null, null);
341 };
342
343 auto set_fieldsplit_preconditioner = [&](auto solver) {
345 SNES snes;
346 CHKERR TSGetSNES(solver, &snes);
347 KSP ksp;
348 CHKERR SNESGetKSP(snes, &ksp);
349 PC pc;
350 CHKERR KSPGetPC(ksp, &pc);
351 PetscBool is_pcfs = PETSC_FALSE;
352 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
353
354 if (is_pcfs == PETSC_TRUE) {
355 auto bc_mng = mField.getInterface<BcManager>();
356 auto name_prb = simple->getProblemName();
357 auto is_all_bc = bc_mng->getBlockIS(name_prb, "ESSENTIAL", "U", 0, 0);
358 int is_all_bc_size;
359 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
360 MOFEM_LOG("EXAMPLE", Sev::inform)
361 << "Field split block size " << is_all_bc_size;
362 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
363 is_all_bc); // boundary block
364 }
366 };
367
368 auto dm = simple->getDM();
369 auto D = smartCreateDMVector(dm);
370 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
371
372 auto solver = pipeline_mng->createTSIM();
373 CHKERR TSSetSolution(solver, D);
374 CHKERR set_time_monitor(dm, solver);
375 CHKERR TSSetSolution(solver, D);
376 CHKERR TSSetFromOptions(solver);
377 CHKERR set_fieldsplit_preconditioner(solver);
378 CHKERR TSSetUp(solver);
379 CHKERR TSSolve(solver, NULL);
380
381 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
382 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
383 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
384
386}
387
390
391 // Processes to set output results are integrated in solveSystem()
392
394}
395
398
407
409}
410
411int main(int argc, char *argv[]) {
412
413 // Initialisation of MoFEM/PETSc and MOAB data structures
414 const char param_file[] = "param_file.petsc";
416
417 // Add logging channel for example
418 auto core_log = logging::core::get();
419 core_log->add_sink(
420 LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
421 LogManager::setLog("EXAMPLE");
422 MOFEM_LOG_TAG("EXAMPLE", "example")
423
424 // Error handling
425 try {
426 // Register MoFEM discrete manager in PETSc
427 DMType dm_name = "DMMOFEM";
428 CHKERR DMRegister_MoFEM(dm_name);
429
430 // Create MOAB instance
431 moab::Core mb_instance; // mesh database
432 moab::Interface &moab = mb_instance; // mesh database interface
433
434 // Create MoFEM instance
435 MoFEM::Core core(moab); // finite element database
436 MoFEM::Interface &m_field = core; // finite element interface
437
438 // Run the main analysis
439 HeatEquation heat_problem(m_field);
440 CHKERR heat_problem.runProgram();
441 }
443
444 // Finish work: cleaning memory, getting statistics, etc.
446
447 return 0;
448}
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEM::FaceElementForcesAndSourcesCore FaceEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
auto integration_rule
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:482
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:545
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
FTensor::Index< 'i', SPACE_DIM > i
constexpr double init_u
int main(int argc, char *argv[])
EntitiesFieldData::EntData EntData
[Define dimension]
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBoundaryTimeScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
static char help[]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpDomainTimesScalarField
constexpr double k
constexpr int SPACE_DIM
[Define dimension]
constexpr double c
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
Definition: helmholtz.cpp:37
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:43
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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: DMMMoFEM.cpp:1015
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
const double D
diffusivity
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:93
constexpr double t
plate stiffness
Definition: plate.cpp:76
static constexpr int approx_order
MoFEM::Interface & mField
MoFEMErrorCode solveSystem()
MoFEMErrorCode outputResults()
HeatEquation(MoFEM::Interface &m_field)
MoFEMErrorCode setupProblem()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode boundaryCondition()
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode runProgram()
MoFEMErrorCode initialCondition()
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
structure for User Loop Methods on finite elements
Basic algebra on fields.
Definition: FieldBlas.hpp:31
Interface for managing meshsets containing materials and boundary conditions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Simple interface for fast problem set-up.
Definition: Simple.hpp:35
PetscInt ts_step
time step
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Monitor solution.
Monitor(SmartPetscObj< DM > dm, boost::shared_ptr< PostProcEle > post_proc)
static constexpr int saveEveryNthStep
SmartPetscObj< DM > dM
MoFEMErrorCode postProcess()
function is run at the end of loop
boost::shared_ptr< PostProcEle > postProc
MoFEMErrorCode operator()()
function is run for every finite element
MoFEMErrorCode preProcess()
function is run at the beginning of loop
Postprocess on face.