v0.13.2
Loading...
Searching...
No Matches
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
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
73 MoFEMErrorCode doWork(int side, EntityType type,
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(
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(
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
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
int main()
Definition: adol-c_atom.cpp:46
constexpr double a
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: DMMoFEM.cpp:491
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:978
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.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:364
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:332
double dt
Definition: heat_method.cpp:26
static char help[]
Definition: heat_method.cpp:24
constexpr int SPACE_DIM
Definition: heat_method.cpp:34
double D
const double v
phase velocity of light in medium (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM > OpGradGrad
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double rho
Definition: plastic.cpp:98
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
boost::shared_ptr< MatrixDouble > uGradPtr
Definition: heat_method.cpp:77
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpRhs(boost::shared_ptr< MatrixDouble > u_grad_ptr)
[Example]
Definition: plastic.cpp:138
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:145
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.
static UId getUniqueIdCalculate(const DofIdx dof, UId ent_uid)
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.
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:279
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:323
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()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
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:27
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:264
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:668
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:476
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.