v0.13.1
heat_method.cpp
Go to the documentation of this file.
1/**
2 * \file heat_method.cpp \example heat_method.cpp
3 *
4 * Calculating geodetic distances using heat method. For details see,
5 *
6 * K. Crane, C. Weischedel, M. Wardetzky, Geodesics in heat: A new approach to
7 * computing distance based on heat flow, ACM Transactions on Graphics , vol.
8 * 32, no. 5, pp. 152:1-152:11, 2013.
9 *
10 * Interent resources:
11 * http://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/
12 * http://www.numerical-tours.com/matlab/meshproc_7_geodesic_poisson/
13 *
14 *
15 * Solving two problems in sequence. Show hot to use form integrators and how to
16 * implement user data operator.
17 *
18 */
19
20#include <MoFEM.hpp>
21
22using namespace MoFEM;
23
24static char help[] = "...\n\n";
25
26double dt = 2; // relative to edge length
27
28#include <BasicFiniteElements.hpp>
29
33
34constexpr int SPACE_DIM = 3;
35
36// Use forms iterators for Grad-Grad term
39
40// Use forms for Mass term
43
44struct Example {
45
46 Example(MoFEM::Interface &m_field) : mField(m_field) {}
47
49
50private:
53
55
65
66 /**
67 * Use problem specific implementation for second stage of heat methid
68 */
69 struct OpRhs : public DomainEleOp {
70
71 OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr);
72
75
76 protected:
77 boost::shared_ptr<MatrixDouble> uGradPtr;
78 };
79};
80
81//! [Run programme]
94}
95//! [Run programme]
96
97//! [Read mesh]
100
104
105 // FIXME: THis part of algorithm is not working in parallel. If you have bit
106 // of free time, feel free to generalise code below.
107
108 Range nodes;
109 CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes);
110 // pinchNodes could be get from BLOCKSET
111 pinchNodes.insert(nodes[0]);
112
113 Range edges;
114 CHKERR mField.get_moab().get_adjacencies(pinchNodes, 1, false, edges,
115 moab::Interface::UNION);
116 double l2;
117 for (auto e : edges) {
118 Range nodes;
119 CHKERR mField.get_moab().get_connectivity(Range(e, e), nodes, false);
120 MatrixDouble coords(nodes.size(), 3);
121 CHKERR mField.get_moab().get_coords(nodes, &coords(0, 0));
122 double l2e = 0;
123 for (int j = 0; j != 3; ++j) {
124 l2e += pow(coords(0, j) - coords(1, j), 2);
125 }
126 l2 = std::max(l2, l2e);
127 }
128
129 dt *= std::sqrt(l2);
130
132}
133//! [Read mesh]
134
135//! [Set up problem]
138
139 // Only one field
141 constexpr int order = 1;
144
146}
147//! [Set up problem]
148
149//! [Set integration rule]
152
153 // Set integration order. To 2p is enough to integrate mass matrix exactly.
154 auto rule = [](int, int, int p) -> int { return 2 * p; };
155
156 // Set integration rule for elements assembling matrix and vector. Note that
157 // rule for vector could be 2p-1, but that not change computation time
158 // significantly. So shorter code is better.
160 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
161 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
162
164}
165//! [Set integration rule]
166
167//! [Create common data]
171}
172//! [Create common data]
173
174//! [Boundary condition]
178}
179//! [Boundary condition]
180
181//! [Push operators to pipeline]
185
186 // This will store gradients at integration points on element
187 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
188
189 // Push element from reference configuration to current configuration in 3d
190 // space
191 auto set_domain_general = [&](auto &pipeline) {
192 auto det_ptr = boost::make_shared<VectorDouble>();
193 auto jac_ptr = boost::make_shared<MatrixDouble>();
194 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
195 pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
196 pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
197 pipeline.push_back(new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
198 };
199
200 // Operators for assembling matrix for first stage
201 auto set_domain_lhs_first = [&](auto &pipeline) {
202 auto one = [](double, double, double) -> double { return 1.; };
203 pipeline.push_back(new OpGradGrad("U", "U", one));
204 auto rho = [](double, double, double) -> double { return 1. / dt; };
205 pipeline.push_back(new OpMass("U", "U", rho));
206 };
207
208 // Nothing is assembled for vector for first stage of heat method. Later
209 // simply value of one is set to elements of right hand side vector at pinch
210 // nodes.
211 auto set_domain_rhs_first = [&](auto &pipeline) {};
212
213 // Operators for assembly of left hand side. This time only Grad-Grad
214 // operator.
215 auto set_domain_lhs_second = [&](auto &pipeline) {
216 auto one = [](double, double, double) { return 1.; };
217 pipeline.push_back(new OpGradGrad("U", "U", one));
218 };
219
220 // Now apply on the right hand side vector X = Grad/||Grad||
221 auto set_domain_rhs_second = [&](auto &pipeline) {
222 pipeline.push_back(new OpCalculateScalarFieldGradient<3>("U", grad_u_ptr));
223 pipeline.push_back(new OpRhs(grad_u_ptr));
224 };
225
226 // Solver first problem
227 auto solve_first = [&]() {
230 auto pipeline_mng = mField.getInterface<PipelineManager>();
231
232 auto solver = pipeline_mng->createKSP();
233 CHKERR KSPSetFromOptions(solver);
234 CHKERR KSPSetUp(solver);
235
236 auto dm = simpleInterface->getDM();
237 auto D = smartCreateDMVector(dm);
238 auto F = smartVectorDuplicate(D);
239
240 // Note add one at nodes which are pinch nodes
241 CHKERR VecZeroEntries(F);
242 auto problem_ptr = mField.get_problem(simple->getProblemName());
243 auto bit_number = mField.get_field_bit_number("U");
244 auto dofs_ptr = problem_ptr->getNumeredRowDofsPtr();
245 for (auto v : pinchNodes) {
246 const auto uid = DofEntity::getUniqueIdCalculate(
247 0, FieldEntity::getLocalUniqueIdCalculate(bit_number, v));
248 auto dof = dofs_ptr->get<Unique_mi_tag>().find(uid);
249 if (dof != dofs_ptr->end())
250 VecSetValue(F, (*dof)->getPetscGlobalDofIdx(), 1, INSERT_VALUES);
251 }
252 CHKERR VecAssemblyBegin(F);
253 CHKERR VecAssemblyEnd(F);
254
255 // Solve problem
256 CHKERR KSPSolve(solver, F, D);
257 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
258 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
259 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
261 };
262
263 // Solve second problem
264 auto solve_second = [&]() {
267
268 // Note that DOFs on pinch nodes are removed from the problem
269 auto prb_mng = mField.getInterface<ProblemsManager>();
270 CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "U",
271 pinchNodes, 0, 1);
272
273 auto *pipeline_mng = mField.getInterface<PipelineManager>();
274
275 auto solver = pipeline_mng->createKSP();
276 CHKERR KSPSetFromOptions(solver);
277 CHKERR KSPSetUp(solver);
278
279 auto dm = simpleInterface->getDM();
280 auto D = smartCreateDMVector(dm);
281 auto F = smartVectorDuplicate(D);
282
283 CHKERR KSPSolve(solver, F, D);
284 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
285 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
286 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
288 };
289
290 // Post-process results
291 auto post_proc = [&](const std::string out_name) {
294 auto tmp_lhs_fe = pipeline_mng->getDomainLhsFE();
295 auto tmp_rhs_fe = pipeline_mng->getDomainRhsFE();
296 auto det_ptr = boost::make_shared<VectorDouble>();
297 auto jac_ptr = boost::make_shared<MatrixDouble>();
298 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
299 pipeline_mng->getDomainLhsFE().reset();
300 pipeline_mng->getDomainRhsFE().reset();
301 auto post_proc_fe =
302 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
303
304 post_proc_fe->getOpPtrVector().push_back(
306 post_proc_fe->getOpPtrVector().push_back(
307 new OpInvertMatrix<SPACE_DIM>(jac_ptr, det_ptr, inv_jac_ptr));
308 post_proc_fe->getOpPtrVector().push_back(
310
311 auto u_ptr = boost::make_shared<VectorDouble>();
312 auto grad_ptr = boost::make_shared<MatrixDouble>();
313
314 post_proc_fe->getOpPtrVector().push_back(
315 new OpCalculateScalarFieldValues("U", u_ptr));
316 post_proc_fe->getOpPtrVector().push_back(
318
320
321 post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
322 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
323
324 {{"U", u_ptr}},
325
326 {{"GRAD_U", grad_ptr}},
327
328 {}, {}));
329
330 pipeline_mng->getDomainRhsFE() = post_proc_fe;
331 CHKERR pipeline_mng->loopFiniteElements();
332 CHKERR post_proc_fe->writeFile(out_name);
333 tmp_lhs_fe = pipeline_mng->getDomainLhsFE() = tmp_lhs_fe;
334 tmp_rhs_fe = pipeline_mng->getDomainRhsFE() = tmp_rhs_fe;
336 };
337
338 set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
339 set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
340 set_domain_lhs_first(pipeline_mng->getOpDomainLhsPipeline());
341 set_domain_rhs_first(pipeline_mng->getOpDomainRhsPipeline());
342
343 CHKERR solve_first();
344 CHKERR post_proc("out_heat_method_first.h5m");
345
346 pipeline_mng->getOpDomainLhsPipeline().clear();
347 pipeline_mng->getOpDomainRhsPipeline().clear();
348
349 set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
350 set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
351 set_domain_lhs_second(pipeline_mng->getOpDomainLhsPipeline());
352 set_domain_rhs_second(pipeline_mng->getOpDomainRhsPipeline());
353
354 CHKERR solve_second();
355 CHKERR post_proc("out_heat_method_second.h5m");
356
358};
359//! [Push operators to pipeline]
360
361//! [Solve]
365}
366//! [Solve]
367
368//! [Postprocess results]
372}
373//! [Postprocess results]
374
375//! [Check results]
379}
380//! [Check results]
381
382int main(int argc, char *argv[]) {
383
384 // Initialisation of MoFEM/PETSc and MOAB data structures
385 const char param_file[] = "param_file.petsc";
387
388 // Add logging channel for example
389 auto core_log = logging::core::get();
390 core_log->add_sink(
391 LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
392 LogManager::setLog("EXAMPLE");
393 MOFEM_LOG_TAG("EXAMPLE", "example");
394
395 try {
396
397 //! [Register MoFEM discrete manager in PETSc]
398 DMType dm_name = "DMMOFEM";
399 CHKERR DMRegister_MoFEM(dm_name);
400 //! [Register MoFEM discrete manager in PETSc
401
402 //! [Create MoAB]
403 moab::Core mb_instance; ///< mesh database
404 moab::Interface &moab = mb_instance; ///< mesh database interface
405 //! [Create MoAB]
406
407 //! [Create MoFEM]
408 MoFEM::Core core(moab); ///< finite element database
409 MoFEM::Interface &m_field = core; ///< finite element database insterface
410 //! [Create MoFEM]
411
412 //! [Example]
413 Example ex(m_field);
414 CHKERR ex.runProblem();
415 //! [Example]
416 }
418
420}
421
422Example::OpRhs::OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr)
423 : DomainEleOp("U", DomainEleOp::OPROW), uGradPtr(u_grad_ptr) {}
424
428
430
431 auto nb_dofs = data.getIndices().size();
432 if (nb_dofs) {
433
434 auto t_grad = getFTensor1FromMat<3>(*uGradPtr);
435
436 auto nb_base_functions = data.getN().size2();
437 auto nb_gauss_pts = getGaussPts().size2();
438 std::array<double, MAX_DOFS_ON_ENTITY> nf;
439 std::fill(nf.begin(), &nf[nb_dofs], 0);
440
441 auto t_diff_base = data.getFTensor1DiffN<3>();
442 auto t_w = getFTensor0IntegrationWeight();
443 auto a = getMeasure();
444
445 for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
446 double alpha = t_w * a;
447
448 const auto l2 = t_grad(i) * t_grad(i);
450 if (l2 > std::numeric_limits<double>::epsilon())
451 t_one(i) = t_grad(i) / std::sqrt(l2);
452 else
453 t_one(i) = 0;
454
455 size_t bb = 0;
456 for (; bb != nb_dofs; ++bb) {
457 nf[bb] -= alpha * t_diff_base(i) * t_one(i);
458 ++t_diff_base;
459 }
460
461 for (; bb < nb_base_functions; ++bb) {
462 ++t_diff_base;
463 }
464
465 ++t_grad;
466 }
467
468 CHKERR VecSetValues<MoFEM::EssentialBcStorage>(getKSPf(), data, &nf[0],
469 ADD_VALUES);
470 }
471
473}
static Index< 'p', 3 > p
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
constexpr double a
MoFEM::FaceElementForcesAndSourcesCore FaceEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
@ 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
#define CHKERR
Inline error check.
Definition: definitions.h:535
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
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:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
FTensor::Index< 'i', SPACE_DIM > i
double dt
Definition: heat_method.cpp:26
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
Definition: heat_method.cpp:32
static char help[]
Definition: heat_method.cpp:24
constexpr int SPACE_DIM
Definition: heat_method.cpp:34
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMass
Definition: heat_method.cpp:42
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 3 > OpGradGrad
Definition: heat_method.cpp:38
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double rho
Definition: plastic.cpp:98
boost::shared_ptr< MatrixDouble > uGradPtr
Definition: heat_method.cpp:77
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpRhs(boost::shared_ptr< MatrixDouble > u_grad_ptr)
[Example]
Definition: plastic.cpp:120
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
Definition: heat_method.cpp:46
Simple * simpleInterface
Definition: helmholtz.cpp:45
Range pinchNodes
Definition: heat_method.cpp:54
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:127
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.