v0.13.0
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 /* This file is part of MoFEM.
21  * MoFEM is free software: you can redistribute it and/or modify it under
22  * the terms of the GNU Lesser General Public License as published by the
23  * Free Software Foundation, either version 3 of the License, or (at your
24  * option) any later version.
25  *
26  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
27  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
28  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
29  * License for more details.
30  *
31  * You should have received a copy of the GNU Lesser General Public
32  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
33 
34 #include <MoFEM.hpp>
35 
36 using namespace MoFEM;
37 
38 static char help[] = "...\n\n";
39 
40 double dt = 2; // relative to edge length
41 
42 #include <BasicFiniteElements.hpp>
43 
47 
48 // Use forms iterators for Grad-Grad term
51 
52 // Use forms for Mass term
55 
56 struct Example {
57 
58  Example(MoFEM::Interface &m_field) : mField(m_field) {}
59 
61 
62 private:
63  MoFEM::Interface &mField;
64  Simple *simpleInterface;
65 
66  Range pinchNodes;
67 
77 
78  /**
79  * Use problem specific implementation for second stage of heat methid
80  */
81  struct OpRhs : public DomainEleOp {
82 
83  OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr);
84 
85  MoFEMErrorCode doWork(int side, EntityType type,
87 
88  protected:
89  boost::shared_ptr<MatrixDouble> uGradPtr;
90  };
91 };
92 
93 //! [Run programme]
96  CHKERR readMesh();
97  CHKERR setupProblem();
98  CHKERR createCommonData();
99  CHKERR boundaryCondition();
100  CHKERR assembleSystem();
101  CHKERR setIntegrationRules();
102  CHKERR solveSystem();
103  CHKERR outputResults();
104  CHKERR checkResults();
106 }
107 //! [Run programme]
108 
109 //! [Read mesh]
112 
113  CHKERR mField.getInterface(simpleInterface);
114  CHKERR simpleInterface->getOptions();
115  CHKERR simpleInterface->loadFile();
116 
117  // FIXME: THis part of algorithm is not working in parallel. If you have bit
118  // of free time, feel free to generalise code below.
119 
120  Range nodes;
121  CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes);
122  // pinchNodes could be get from BLOCKSET
123  pinchNodes.insert(nodes[0]);
124 
125  Range edges;
126  CHKERR mField.get_moab().get_adjacencies(pinchNodes, 1, false, edges,
127  moab::Interface::UNION);
128  double l2;
129  for (auto e : edges) {
130  Range nodes;
131  CHKERR mField.get_moab().get_connectivity(Range(e, e), nodes, false);
132  MatrixDouble coords(nodes.size(), 3);
133  CHKERR mField.get_moab().get_coords(nodes, &coords(0, 0));
134  double l2e = 0;
135  for (int j = 0; j != 3; ++j) {
136  l2e += pow(coords(0, j) - coords(1, j), 2);
137  }
138  l2 = std::max(l2, l2e);
139  }
140 
141  dt *= std::sqrt(l2);
142 
144 }
145 //! [Read mesh]
146 
147 //! [Set up problem]
150 
151  // Only one field
152  CHKERR simpleInterface->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
153  constexpr int order = 1;
154  CHKERR simpleInterface->setFieldOrder("U", order);
155  CHKERR simpleInterface->setUp();
156 
158 }
159 //! [Set up problem]
160 
161 //! [Set integration rule]
164 
165  // Set integration order. To 2p is enough to integrate mass matrix exactly.
166  auto rule = [](int, int, int p) -> int { return 2 * p; };
167 
168  // Set integration rule for elements assembling matrix and vector. Note that
169  // rule for vector could be 2p-1, but that not change computation time
170  // significantly. So shorter code is better.
171  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
172  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
173  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
174 
176 }
177 //! [Set integration rule]
178 
179 //! [Create common data]
183 }
184 //! [Create common data]
185 
186 //! [Boundary condition]
190 }
191 //! [Boundary condition]
192 
193 //! [Push operators to pipeline]
196  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
197 
198  // This will store gradients at integration points on element
199  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
200 
201  // Push element from reference configuration to current configuration in 3d
202  // space
203  auto set_domain_general = [&](auto &pipeline) {
204  auto det_ptr = boost::make_shared<VectorDouble>();
205  auto jac_ptr = boost::make_shared<MatrixDouble>();
206  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
207  pipeline.push_back(new OpCalculateHOJacForFaceEmbeddedIn3DSpace(jac_ptr));
208  pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
209  pipeline.push_back(new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
210  };
211 
212  // Operators for assembling matrix for first stage
213  auto set_domain_lhs_first = [&](auto &pipeline) {
214  auto one = [](double, double, double) -> double { return 1.; };
215  pipeline.push_back(new OpGradGrad("U", "U", one));
216  auto rho = [](double, double, double) -> double { return 1. / dt; };
217  pipeline.push_back(new OpMass("U", "U", rho));
218  };
219 
220  // Nothing is assembled for vector for first stage of heat method. Later
221  // simply value of one is set to elements of right hand side vector at pinch
222  // nodes.
223  auto set_domain_rhs_first = [&](auto &pipeline) {};
224 
225  // Operators for assembly of left hand side. This time only Grad-Grad
226  // operator.
227  auto set_domain_lhs_second = [&](auto &pipeline) {
228  auto one = [](double, double, double) { return 1.; };
229  pipeline.push_back(new OpGradGrad("U", "U", one));
230  };
231 
232  // Now apply on the right hand side vector X = Grad/||Grad||
233  auto set_domain_rhs_second = [&](auto &pipeline) {
234  pipeline.push_back(new OpCalculateScalarFieldGradient<3>("U", grad_u_ptr));
235  pipeline.push_back(new OpRhs(grad_u_ptr));
236  };
237 
238  // Solver first problem
239  auto solve_first = [&]() {
241  auto simple = mField.getInterface<Simple>();
242  auto pipeline_mng = mField.getInterface<PipelineManager>();
243 
244  auto solver = pipeline_mng->createKSP();
245  CHKERR KSPSetFromOptions(solver);
246  CHKERR KSPSetUp(solver);
247 
248  auto dm = simpleInterface->getDM();
249  auto D = smartCreateDMVector(dm);
250  auto F = smartVectorDuplicate(D);
251 
252  // Note add one at nodes which are pinch nodes
253  CHKERR VecZeroEntries(F);
254  auto problem_ptr = mField.get_problem(simple->getProblemName());
255  auto bit_number = mField.get_field_bit_number("U");
256  auto dofs_ptr = problem_ptr->getNumeredRowDofsPtr();
257  for (auto v : pinchNodes) {
258  const auto uid = DofEntity::getUniqueIdCalculate(
259  0, FieldEntity::getLocalUniqueIdCalculate(bit_number, v));
260  auto dof = dofs_ptr->get<Unique_mi_tag>().find(uid);
261  if (dof != dofs_ptr->end())
262  VecSetValue(F, (*dof)->getPetscGlobalDofIdx(), 1, INSERT_VALUES);
263  }
264  CHKERR VecAssemblyBegin(F);
265  CHKERR VecAssemblyEnd(F);
266 
267  // Solve problem
268  CHKERR KSPSolve(solver, F, D);
269  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
270  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
271  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
273  };
274 
275  // Solve second problem
276  auto solve_second = [&]() {
278  auto simple = mField.getInterface<Simple>();
279 
280  // Note that DOFs on pinch nodes are removed from the problem
281  auto prb_mng = mField.getInterface<ProblemsManager>();
282  CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(), "U",
283  pinchNodes, 0, 1);
284 
285  auto *pipeline_mng = mField.getInterface<PipelineManager>();
286 
287  auto solver = pipeline_mng->createKSP();
288  CHKERR KSPSetFromOptions(solver);
289  CHKERR KSPSetUp(solver);
290 
291  auto dm = simpleInterface->getDM();
292  auto D = smartCreateDMVector(dm);
293  auto F = smartVectorDuplicate(D);
294 
295  CHKERR KSPSolve(solver, F, D);
296  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
297  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
298  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
300  };
301 
302  // Post-process results
303  auto post_proc = [&](const std::string out_name) {
304  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
306  auto tmp_lhs_fe = pipeline_mng->getDomainLhsFE();
307  auto tmp_rhs_fe = pipeline_mng->getDomainRhsFE();
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  pipeline_mng->getDomainLhsFE().reset();
312  pipeline_mng->getDomainRhsFE().reset();
313  auto post_proc_fe =
314  boost::make_shared<PostProcFaceOnRefinedMeshFor2D>(mField);
315  post_proc_fe->generateReferenceElementMesh();
316  post_proc_fe->getOpPtrVector().push_back(
318  post_proc_fe->getOpPtrVector().push_back(
319  new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
320  post_proc_fe->getOpPtrVector().push_back(
321  new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
322  post_proc_fe->addFieldValuesPostProc("U");
323  post_proc_fe->addFieldValuesGradientPostProc("U");
324  pipeline_mng->getDomainRhsFE() = post_proc_fe;
325  CHKERR pipeline_mng->loopFiniteElements();
326  CHKERR post_proc_fe->writeFile(out_name);
327  tmp_lhs_fe = pipeline_mng->getDomainLhsFE() = tmp_lhs_fe;
328  tmp_rhs_fe = pipeline_mng->getDomainRhsFE() = tmp_rhs_fe;
330  };
331 
332  set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
333  set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
334  set_domain_lhs_first(pipeline_mng->getOpDomainLhsPipeline());
335  set_domain_rhs_first(pipeline_mng->getOpDomainRhsPipeline());
336 
337  CHKERR solve_first();
338  CHKERR post_proc("out_heat_method_first.h5m");
339 
340  pipeline_mng->getOpDomainLhsPipeline().clear();
341  pipeline_mng->getOpDomainRhsPipeline().clear();
342 
343  set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
344  set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
345  set_domain_lhs_second(pipeline_mng->getOpDomainLhsPipeline());
346  set_domain_rhs_second(pipeline_mng->getOpDomainRhsPipeline());
347 
348  CHKERR solve_second();
349  CHKERR post_proc("out_heat_method_second.h5m");
350 
352 };
353 //! [Push operators to pipeline]
354 
355 //! [Solve]
359 }
360 //! [Solve]
361 
362 //! [Postprocess results]
366 }
367 //! [Postprocess results]
368 
369 //! [Check results]
373 }
374 //! [Check results]
375 
376 int main(int argc, char *argv[]) {
377 
378  // Initialisation of MoFEM/PETSc and MOAB data structures
379  const char param_file[] = "param_file.petsc";
380  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
381 
382  // Add logging channel for example
383  auto core_log = logging::core::get();
384  core_log->add_sink(
386  LogManager::setLog("EXAMPLE");
387  MOFEM_LOG_TAG("EXAMPLE", "example");
388 
389  try {
390 
391  //! [Register MoFEM discrete manager in PETSc]
392  DMType dm_name = "DMMOFEM";
393  CHKERR DMRegister_MoFEM(dm_name);
394  //! [Register MoFEM discrete manager in PETSc
395 
396  //! [Create MoAB]
397  moab::Core mb_instance; ///< mesh database
398  moab::Interface &moab = mb_instance; ///< mesh database interface
399  //! [Create MoAB]
400 
401  //! [Create MoFEM]
402  MoFEM::Core core(moab); ///< finite element database
403  MoFEM::Interface &m_field = core; ///< finite element database insterface
404  //! [Create MoFEM]
405 
406  //! [Example]
407  Example ex(m_field);
408  CHKERR ex.runProblem();
409  //! [Example]
410  }
411  CATCH_ERRORS;
412 
414 }
415 
416 Example::OpRhs::OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr)
417  : DomainEleOp("U", DomainEleOp::OPROW), uGradPtr(u_grad_ptr) {}
418 
422 
424 
425  auto nb_dofs = data.getIndices().size();
426  if (nb_dofs) {
427 
428  auto t_grad = getFTensor1FromMat<3>(*uGradPtr);
429 
430  auto nb_base_functions = data.getN().size2();
431  auto nb_gauss_pts = getGaussPts().size2();
432  std::array<double, MAX_DOFS_ON_ENTITY> nf;
433  std::fill(nf.begin(), &nf[nb_dofs], 0);
434 
435  auto t_diff_base = data.getFTensor1DiffN<3>();
436  auto t_w = getFTensor0IntegrationWeight();
437  auto a = getMeasure();
438 
439  for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
440  double alpha = t_w * a;
441 
442  const auto l2 = t_grad(i) * t_grad(i);
444  if (l2 > std::numeric_limits<double>::epsilon())
445  t_one(i) = t_grad(i) / std::sqrt(l2);
446  else
447  t_one(i) = 0;
448 
449  size_t bb = 0;
450  for (; bb != nb_dofs; ++bb) {
451  nf[bb] -= alpha * t_diff_base(i) * t_one(i);
452  ++t_diff_base;
453  }
454 
455  for (; bb < nb_base_functions; ++bb) {
456  ++t_diff_base;
457  }
458 
459  ++t_grad;
460  }
461 
462  CHKERR VecSetValues<MoFEM::EssentialBcStorage>(getKSPf(), data, &nf[0],
463  ADD_VALUES);
464  }
465 
467 }
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
constexpr double a
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
constexpr double rho
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
constexpr double alpha
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ 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
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
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:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
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:342
FTensor::Index< 'i', SPACE_DIM > i
double dt
Definition: heat_method.cpp:40
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
Definition: heat_method.cpp:46
static char help[]
Definition: heat_method.cpp:38
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMass
Definition: heat_method.cpp:54
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 3 > OpGradGrad
Definition: heat_method.cpp:50
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:88
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1096
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
boost::shared_ptr< MatrixDouble > uGradPtr
Definition: heat_method.cpp:89
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]
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
Definition: heat_method.cpp:58
Range pinchNodes
Definition: heat_method.cpp:66
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
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.
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.
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:33
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator