v0.13.1
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 105 of file heat_equation.cpp.

Constructor & Destructor Documentation

◆ HeatEquation()

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

Definition at line 130 of file heat_equation.cpp.

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

Member Function Documentation

◆ assembleSystem()

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

Definition at line 219 of file heat_equation.cpp.

219 {
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}
@ 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
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 1 > OpBoundaryTimeScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalarField< 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:76
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
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()
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.

◆ boundaryCondition()

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

Definition at line 197 of file heat_equation.cpp.

197 {
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}
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:548
FTensor::Index< 'i', SPACE_DIM > i
Simple interface for fast problem set-up.
Definition: BcManager.hpp:27
Simple interface for fast problem set-up.
Definition: Simple.hpp:35

◆ initialCondition()

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

Definition at line 175 of file heat_equation.cpp.

175 {
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}
@ BLOCKSET
Definition: definitions.h:161
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:31
Interface for managing meshsets containing materials and boundary conditions.

◆ outputResults()

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

Definition at line 388 of file heat_equation.cpp.

388 {
390
391 // Processes to set output results are integrated in solveSystem()
392
394}

◆ readMesh()

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

Definition at line 132 of file heat_equation.cpp.

132 {
134
135 auto *simple = mField.getInterface<Simple>();
137 CHKERR simple->getOptions();
138 CHKERR simple->loadFile();
139
141}

◆ runProgram()

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

Definition at line 396 of file heat_equation.cpp.

396 {
398
407
409}
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 159 of file heat_equation.cpp.

159 {
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}
auto integration_rule
static constexpr int approx_order

◆ setupProblem()

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

Definition at line 143 of file heat_equation.cpp.

143 {
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}
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)

◆ solveSystem()

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

Definition at line 312 of file heat_equation.cpp.

312 {
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}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
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
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:977
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
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
const double D
diffusivity
Monitor solution.

Member Data Documentation

◆ boundaryMarker

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

Definition at line 127 of file heat_equation.cpp.

◆ mField

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

Definition at line 124 of file heat_equation.cpp.


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