v0.14.0
Loading...
Searching...
No Matches
Public Member Functions | Private Member Functions | Private Attributes | List of all members
HeatEquation Struct Reference
Collaboration diagram for HeatEquation:
[legend]

Public Member Functions

 HeatEquation (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProgram ()
 

Private Member Functions

MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 
MoFEMErrorCode initialCondition ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 

Private Attributes

MoFEM::InterfacemField
 
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
 

Detailed Description

Examples
heat_equation.cpp.

Definition at line 91 of file heat_equation.cpp.

Constructor & Destructor Documentation

◆ HeatEquation()

HeatEquation::HeatEquation ( MoFEM::Interface m_field)
Examples
heat_equation.cpp.

Definition at line 116 of file heat_equation.cpp.

116: mField(m_field) {}
MoFEM::Interface & mField

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode HeatEquation::assembleSystem ( )
private
Examples
heat_equation.cpp.

Definition at line 205 of file heat_equation.cpp.

205 {
207
208 auto add_domain_lhs_ops = [&](auto &pipeline) {
210 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
211 pipeline.push_back(new OpDomainGradGrad(
212 "U", "U", [](double, double, double) -> double { return k; }));
213 auto get_c = [this](const double, const double, const double) {
214 auto pipeline_mng = mField.getInterface<PipelineManager>();
215 auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
216 return c * fe_domain_lhs->ts_a;
217 };
218 pipeline.push_back(new OpDomainMass("U", "U", get_c));
219 pipeline.push_back(new OpUnSetBc("U"));
220 };
221
222 auto add_domain_rhs_ops = [&](auto &pipeline) {
224 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
225 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
226 auto dot_u_at_gauss_pts = boost::make_shared<VectorDouble>();
227 pipeline.push_back(new OpCalculateScalarFieldGradient<SPACE_DIM>(
228 "U", grad_u_at_gauss_pts));
229 pipeline.push_back(
230 new OpCalculateScalarFieldValuesDot("U", dot_u_at_gauss_pts));
231 pipeline.push_back(new OpDomainGradTimesVec(
232 "U", grad_u_at_gauss_pts,
233 [](double, double, double) -> double { return k; }));
234 pipeline.push_back(new OpDomainTimesScalarField(
235 "U", dot_u_at_gauss_pts,
236 [](const double, const double, const double) { return c; }));
237 auto source_term = [&](const double x, const double y, const double z) {
238 auto pipeline_mng = mField.getInterface<PipelineManager>();
239 auto &fe_domain_lhs = pipeline_mng->getDomainRhsFE();
240 const auto t = fe_domain_lhs->ts_t;
241 return 1e1 * pow(M_E, -M_PI * M_PI * t) * sin(1. * M_PI * x) *
242 sin(2. * M_PI * y);
243 };
244 pipeline.push_back(new OpDomainSource("U", source_term));
245 pipeline.push_back(new OpUnSetBc("U"));
246 };
247
248 auto add_boundary_lhs_ops = [&](auto &pipeline) {
249 pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
250 pipeline.push_back(new OpBoundaryMass(
251 "U", "U", [](const double, const double, const double) { return c; }));
252 pipeline.push_back(new OpUnSetBc("U"));
253 };
254 auto add_boundary_rhs_ops = [&](auto &pipeline) {
255 pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
256 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
257 auto boundary_function = [&](const double x, const double y,
258 const double z) {
259 auto pipeline_mng = mField.getInterface<PipelineManager>();
260 auto &fe_rhs = pipeline_mng->getBoundaryRhsFE();
261 const auto t = fe_rhs->ts_t;
262 return 0;
263 // abs(0.1 * pow(M_E, -M_PI * M_PI * t) * sin(2. * M_PI * x) *
264 // sin(3. * M_PI * y));
265 };
266 pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
267 pipeline.push_back(new OpBoundaryTimeScalarField(
268 "U", u_at_gauss_pts,
269 [](const double, const double, const double) { return c; }));
270 pipeline.push_back(new OpBoundarySource("U", boundary_function));
271 pipeline.push_back(new OpUnSetBc("U"));
272 };
273
274 auto pipeline_mng = mField.getInterface<PipelineManager>();
275 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
276 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
277 add_boundary_lhs_ops(pipeline_mng->getOpBoundaryLhsPipeline());
278 add_boundary_rhs_ops(pipeline_mng->getOpBoundaryRhsPipeline());
279
281}
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
constexpr double k
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
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_T > OpCalculateScalarFieldValuesDot
constexpr double t
plate stiffness
Definition: plate.cpp:59
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Add operators pushing bases from local to physical configuration.
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.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ boundaryCondition()

MoFEMErrorCode HeatEquation::boundaryCondition ( )
private
Examples
heat_equation.cpp.

Definition at line 183 of file heat_equation.cpp.

183 {
185
186 auto bc_mng = mField.getInterface<BcManager>();
187 auto *simple = mField.getInterface<Simple>();
188 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ESSENTIAL",
189 "U", 0, 0);
190
191 auto &bc_map = bc_mng->getBcMapByBlockName();
192 boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
193 for (auto b : bc_map) {
194 if (std::regex_match(b.first, std::regex("(.*)ESSENTIAL(.*)"))) {
195 boundaryMarker->resize(b.second->bcMarkers.size(), 0);
196 for (int i = 0; i != b.second->bcMarkers.size(); ++i) {
197 (*boundaryMarker)[i] |= b.second->bcMarkers[i];
198 }
199 }
200 }
201
203}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'i', SPACE_DIM > i
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Simple interface for fast problem set-up.
Definition: Simple.hpp:27

◆ initialCondition()

MoFEMErrorCode HeatEquation::initialCondition ( )
private
Examples
heat_equation.cpp.

Definition at line 161 of file heat_equation.cpp.

161 {
163
164 // Get surface entities form blockset, set initial values in those
165 // blocksets. To keep it simple, it is assumed that inital values are on
166 // blockset 1
168 Range inner_surface;
169 CHKERR mField.getInterface<MeshsetsManager>()->getEntitiesByDimension(
170 1, BLOCKSET, 2, inner_surface, true);
171 if (!inner_surface.empty()) {
172 Range inner_surface_verts;
173 CHKERR mField.get_moab().get_connectivity(inner_surface,
174 inner_surface_verts, false);
175 CHKERR mField.getInterface<FieldBlas>()->setField(
176 init_u, MBVERTEX, inner_surface_verts, "U");
177 }
178 }
179
181}
@ BLOCKSET
Definition: definitions.h:148
bool checkMeshset(const int ms_id, const CubitBCType cubit_bc_type) const
check for CUBIT Id and CUBIT type
constexpr double init_u
virtual moab::Interface & get_moab()=0
Basic algebra on fields.
Definition: FieldBlas.hpp:21
Interface for managing meshsets containing materials and boundary conditions.

◆ outputResults()

MoFEMErrorCode HeatEquation::outputResults ( )
private
Examples
heat_equation.cpp.

Definition at line 416 of file heat_equation.cpp.

416 {
418
419 // Processes to set output results are integrated in solveSystem()
420
422}

◆ readMesh()

MoFEMErrorCode HeatEquation::readMesh ( )
private
Examples
heat_equation.cpp.

Definition at line 118 of file heat_equation.cpp.

118 {
120
121 auto *simple = mField.getInterface<Simple>();
122 CHKERR mField.getInterface(simple);
123 CHKERR simple->getOptions();
124 CHKERR simple->loadFile();
125
127}

◆ runProgram()

MoFEMErrorCode HeatEquation::runProgram ( )
Examples
heat_equation.cpp.

Definition at line 424 of file heat_equation.cpp.

424 {
426
435
437}
MoFEMErrorCode solveSystem()
MoFEMErrorCode outputResults()
MoFEMErrorCode setupProblem()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode initialCondition()

◆ setIntegrationRules()

MoFEMErrorCode HeatEquation::setIntegrationRules ( )
private
Examples
heat_equation.cpp.

Definition at line 145 of file heat_equation.cpp.

145 {
147
148 auto integration_rule = [](int o_row, int o_col, int approx_order) {
149 return 2 * approx_order;
150 };
151
152 auto *pipeline_mng = mField.getInterface<PipelineManager>();
153 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
154 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
155 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
156 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
157
159}
auto integration_rule
static constexpr int approx_order

◆ setupProblem()

MoFEMErrorCode HeatEquation::setupProblem ( )
private
Examples
heat_equation.cpp.

Definition at line 129 of file heat_equation.cpp.

129 {
131
132 auto *simple = mField.getInterface<Simple>();
133 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
134 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
135
136 int order = 3;
137 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
138 CHKERR simple->setFieldOrder("U", order);
139
140 CHKERR simple->setUp();
141
143}
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
constexpr int order
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)

◆ solveSystem()

MoFEMErrorCode HeatEquation::solveSystem ( )
private

That to work, you have to create solver, as follows,

auto solver = // pipeline_mng->createTSIM( simple->getDM());

That is explicitly use use Simple DM to create solver for DM. Pipeline menage by default creat copy of DM, in case several solvers are used the same DM.

Alternatively you can get dm directly from the solver, i.e.

DM ts_dm;
CHKERR TSGetDM(solver, &ts_dm);
CHKERR DMTSSetIJacobian(
ts_dm, CalcJacobian::set, getDMTsCtx(ts_dm).get());
static PetscErrorCode set(TS ts, PetscReal t, Vec u, Vec u_t, PetscReal a, Mat A, Mat B, void *ctx)
Examples
heat_equation.cpp.

Definition at line 300 of file heat_equation.cpp.

300 {
302
303 auto *simple = mField.getInterface<Simple>();
304 auto *pipeline_mng = mField.getInterface<PipelineManager>();
305
306 auto create_post_process_element = [&]() {
307 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
308 auto det_ptr = boost::make_shared<VectorDouble>();
309 auto jac_ptr = boost::make_shared<MatrixDouble>();
310 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
311 post_proc_fe->getOpPtrVector().push_back(new OpCalculateHOJac<2>(jac_ptr));
312 post_proc_fe->getOpPtrVector().push_back(
313 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
314 post_proc_fe->getOpPtrVector().push_back(
315 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
316
317 auto u_ptr = boost::make_shared<VectorDouble>();
318 post_proc_fe->getOpPtrVector().push_back(
319 new OpCalculateScalarFieldValues("U", u_ptr));
320
322
323 post_proc_fe->getOpPtrVector().push_back(
324
325 new OpPPMap(post_proc_fe->getPostProcMesh(),
326 post_proc_fe->getMapGaussPts(),
327
328 {{"U", u_ptr}},
329
330 {}, {}, {}
331
332 )
333
334 );
335
336 return post_proc_fe;
337 };
338
339 auto set_time_monitor = [&](auto dm, auto solver) {
341 boost::shared_ptr<Monitor> monitor_ptr(
342 new Monitor(dm, create_post_process_element()));
343 boost::shared_ptr<ForcesAndSourcesCore> null;
344 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
345 monitor_ptr, null, null);
347 };
348
349 auto set_fieldsplit_preconditioner = [&](auto solver) {
351 SNES snes;
352 CHKERR TSGetSNES(solver, &snes);
353 KSP ksp;
354 CHKERR SNESGetKSP(snes, &ksp);
355 PC pc;
356 CHKERR KSPGetPC(ksp, &pc);
357 PetscBool is_pcfs = PETSC_FALSE;
358 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
359
360 if (is_pcfs == PETSC_TRUE) {
361 auto bc_mng = mField.getInterface<BcManager>();
362 auto name_prb = simple->getProblemName();
363 auto is_all_bc = bc_mng->getBlockIS(name_prb, "ESSENTIAL", "U", 0, 0);
364 int is_all_bc_size;
365 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
366 MOFEM_LOG("EXAMPLE", Sev::inform)
367 << "Field split block size " << is_all_bc_size;
368 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
369 is_all_bc); // boundary block
370 }
372 };
373
374 /**
375 * That to work, you have to create solver, as follows,
376 \code
377 auto solver = // pipeline_mng->createTSIM( simple->getDM());
378 \endcode
379 That is explicitly use use Simple DM to create solver for DM. Pipeline
380 menage by default creat copy of DM, in case several solvers are used the
381 same DM.
382
383 Alternatively you can get dm directly from the solver, i.e.
384 \code
385 DM ts_dm;
386 CHKERR TSGetDM(solver, &ts_dm);
387 CHKERR DMTSSetIJacobian(
388 ts_dm, CalcJacobian::set, getDMTsCtx(ts_dm).get());
389 \endcode
390 */
391 auto set_user_ts_jacobian = [&](auto dm) {
393 CHKERR DMTSSetIJacobian(dm, CalcJacobian::set, getDMTsCtx(dm).get());
395 };
396
397 auto dm = simple->getDM();
398 auto D = createDMVector(dm);
399 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
400
401 auto solver = pipeline_mng->createTSIM(
402 simple->getDM()); // Note DM is set as argument. If DM is not, internal
403 // copy of pipeline DM is created.
404 CHKERR set_user_ts_jacobian(dm);
405 CHKERR set_time_monitor(dm, solver);
406 CHKERR TSSetSolution(solver, D);
407 CHKERR TSSetFromOptions(solver);
408 CHKERR set_fieldsplit_preconditioner(solver);
409 CHKERR TSSetUp(solver);
410
411 CHKERR TSSolve(solver, D);
412
414}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
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:509
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
double D
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:1042
auto getDMTsCtx(DM dm)
Get TS context data structure used by DM.
Definition: DMMoFEM.hpp:1045
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Post post-proc data at points from hash maps.
Set inverse jacobian to base functions.
[Push operators to pipeline]

Member Data Documentation

◆ boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > HeatEquation::boundaryMarker
private
Examples
heat_equation.cpp.

Definition at line 113 of file heat_equation.cpp.

◆ mField

MoFEM::Interface& HeatEquation::mField
private
Examples
heat_equation.cpp.

Definition at line 110 of file heat_equation.cpp.


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