v0.14.0
phase.cpp
Go to the documentation of this file.
1 /**
2  * \file phase.cpp
3  * \example phase.cpp
4  *
5  * Calculate light phase using <a
6  * href="https://en.wikipedia.org/wiki/Transport-of-intensity_equation">Transport
7  * intensity equation</a>.
8  *
9  * See paper \cite parvizi2015practical
10  *
11  */
12
13 #include <MoFEM.hpp>
14
15 using namespace MoFEM;
16
17 static char help[] = "...\n\n";
18
19 #include <BasicFiniteElements.hpp>
20
25
26 constexpr std::array<int, 3> d1_savitzky_golay_w3_p2 = {-1, 0, 1};
27 constexpr std::array<int, 5> d1_savitzky_golay_w5_p2 = {-2, -1, 0, 1, 2};
28 constexpr std::array<int, 7> d1_savitzky_golay_w7_p2 = {-3, -2, -1, 0, 1, 2, 3};
29 constexpr std::array<int, 9> d1_savitzky_golay_w9_p2 = {-4, -3, -2, -1, 0,
30  1, 2, 3, 4};
31
32 constexpr std::array<int, 5> d1_savitzky_golay_w5_p4 = {1, -8, 0, 8, -1};
33 constexpr std::array<int, 7> d1_savitzky_golay_w7_p4 = {22, -67, -58, 0,
34  58, 67, -22};
35 constexpr std::array<int, 9> d1_savitzky_golay_w9_p4 = {
36  86, -142, -193, -126, 0, 126, 193, 142, -86};
37
38 constexpr std::array<int, 10> d1_normalisation_p2 = {0, 0, 0, 2, 0,
39  10, 0, 28, 0, 60};
40 constexpr std::array<int, 10> d1_normalisation_p4 = {0, 0, 0, 0, 0,
41  12, 0, 252, 0, 1188};
42
43 const int *d1_savitzky_golay_p2[] = {nullptr, nullptr,
44  nullptr, d1_savitzky_golay_w3_p2.data(),
45  nullptr, d1_savitzky_golay_w5_p2.data(),
46  nullptr, d1_savitzky_golay_w7_p2.data(),
47  nullptr, d1_savitzky_golay_w9_p2.data()};
48
49 const int *d1_savitzky_golay_p4[] = {nullptr,
50  nullptr,
51  nullptr,
52  nullptr,
53  nullptr,
55  nullptr,
57  nullptr,
59
60 const int **d1_savitzky_golay[] = {nullptr, nullptr, d1_savitzky_golay_p2,
61  nullptr, d1_savitzky_golay_p4};
62 const int *d1_normalisation[] = {nullptr, nullptr, d1_normalisation_p2.data(),
63  nullptr, d1_normalisation_p4.data()};
64
71  GAUSS>::OpMixDivTimesScalar<2>;
73  PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
74
75 static double k = 1; // wave number
76 static int order_savitzky_golay = 2;
77 static int window_savitzky_golay = 3;
78
79 struct Example {
80
81  Example(MoFEM::Interface &m_field) : mField(m_field) {}
82
83  MoFEMErrorCode runProblem();
84
85 private:
86  MoFEM::Interface &mField;
87  Simple *simpleInterface;
88  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
89
91  MoFEMErrorCode setupProblem();
92  MoFEMErrorCode assembleSystemIntensity();
93  MoFEMErrorCode assembleSystemFlux();
94  MoFEMErrorCode solveSystem();
95  MoFEMErrorCode calculateFlux(double &calc_flux);
96  MoFEMErrorCode outputResults(const int i);
97
98  enum BoundingBox {
99  CENTER_X = 0,
105  LAST_BB
106  };
107
108  static std::vector<double> rZ;
109  static std::vector<MatrixInt> iI;
110  static std::array<double, LAST_BB> aveMaxMin;
111
112  static int focalIndex;
114  static const int *savitzkyGolayWeights;
115
116  static std::pair<int, int> getCoordsInImage(double x, double y);
117
118  static double rhsSource(const double x, const double y, const double);
119  static double lhsFlux(const double x, const double y, const double);
120
121  struct BoundaryOp;
122 };
123
124 std::vector<double> Example::rZ;
125 std::vector<MatrixInt> Example::iI;
126 std::array<double, Example::BoundingBox::LAST_BB> Example::aveMaxMin;
130
131 std::pair<int, int> Example::getCoordsInImage(double x, double y) {
132
133  auto &m = iI[focalIndex];
134  x -= aveMaxMin[MIN_X];
135  y -= aveMaxMin[MIN_Y];
136  x *= (m.size1() - 1) / (aveMaxMin[MAX_X] - aveMaxMin[MIN_X]);
137  y *= (m.size2() - 1) / (aveMaxMin[MAX_Y] - aveMaxMin[MIN_Y]);
138  const auto p = std::make_pair<int, int>(std::round(x), std::round(y));
139
140 #ifndef NDEBUG
141  if (p.first < 0 && p.first >= m.size1())
142  THROW_MESSAGE("Wrong index");
143  if (p.second < 0 && p.second >= m.size2())
144  THROW_MESSAGE("Wrong index");
145 #endif
146
147  return p;
148 }
149
150 double Example::rhsSource(const double x, const double y, const double) {
151  const auto idx = getCoordsInImage(x, y);
152
153  double v = 0;
154  for (auto w = 0; w != window_savitzky_golay; ++w) {
155  const auto i = focalIndex - (window_savitzky_golay - 1) / 2 + w;
156  const auto &intensity = iI[i];
157  v += intensity(idx.first, idx.second) * savitzkyGolayWeights[w];
158  }
159  v = static_cast<double>(v) / savitzkyGolayNormalisation;
160
161  const auto dz = rZ[focalIndex + 1] - rZ[focalIndex - 1];
162  return -k * v / dz;
163 }
164
165 double Example::lhsFlux(const double x, const double y, const double) {
166  const auto idx = getCoordsInImage(x, y);
167  const auto &m = iI[focalIndex];
168  return 1. / m(idx.first, idx.second);
169 }
170
172  BoundaryOp(boost::shared_ptr<MatrixDouble> flux_ptr, double &glob_flux)
173  : EdgeEleOp(NOSPACE, OPSPACE), fluxPtr(flux_ptr), globFlux(glob_flux) {}
174
175  MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
177  const auto nb_gauss_pts = getGaussPts().size2();
178  auto t_flux = getFTensor1FromMat<3>(*fluxPtr);
179  auto t_normal = getFTensor1Normal();
180  auto t_w = getFTensor0IntegrationWeight();
182  for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
183  globFlux += t_w * t_normal(i) * t_flux(i);
184  ++t_flux;
185  ++t_w;
186  };
188  }
189
190 private:
191  boost::shared_ptr<MatrixDouble> fluxPtr;
192  double &globFlux;
193 };
194
195 //! [run problem]
199  CHKERR setupProblem();
200
201  char images_array_files[255] = "out_arrays.txt";
202  CHKERR PetscOptionsGetString(PETSC_NULL, "", "-images_arrays",
203  images_array_files, 255, PETSC_NULL);
204
205  // Look into "extract_data.ipynb" file, it is used to generate data file with
206  // combined stack of the images.
207  auto read_images = [&]() {
208  std::ifstream in;
209  in.open(images_array_files);
210  std::vector<int> values;
211  values.insert(values.end(), std::istream_iterator<int>(in),
212  std::istream_iterator<int>());
213  MOFEM_LOG("WORLD", Sev::inform) << "Read data size " << values.size();
214  in.close();
215  return values;
216  };
217
218  auto structure_data = [&](auto &&data) {
219  constexpr double scale = 1e4; // scale to float
220  auto it = data.begin();
221  if (it == data.end()) {
222  THROW_MESSAGE("No images");
223  }
224  rZ.reserve(*it);
225  iI.reserve(*it);
226  MOFEM_LOG("WORLD", Sev::inform) << "Number of images " << *it;
227  it++;
228  for (; it != data.end();) {
229  rZ.emplace_back(*(it++) / scale);
230  const auto r = *(it++);
231  const auto c = *(it++);
232  iI.emplace_back(r, c);
233  MOFEM_LOG("WORLD", Sev::inform)
234  << "Read data set " << rZ.back() << " size " << r << " by " << c;
235  auto &m = iI.back();
236  for (auto rit = m.begin1(); rit != m.end1(); ++rit) {
237  for (auto cit = rit.begin(); cit != rit.end(); ++cit) {
238  *cit = *(it++);
239  }
240  }
241  }
242  };
243
245
246  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-k", &k, PETSC_NULL);
247  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order_savitzky_golay",
248  &order_savitzky_golay, PETSC_NULL);
249  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-window_savitzky_golay",
250  &window_savitzky_golay, PETSC_NULL);
251
252  int window_shift = 0;
253  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-window_shift", &window_shift,
254  PETSC_NULL);
255
256  MOFEM_LOG("WORLD", Sev::inform) << "Wave number " << k;
257  MOFEM_LOG("WORLD", Sev::inform)
258  << "Order Savitzky Golay " << order_savitzky_golay;
259  MOFEM_LOG("WORLD", Sev::inform)
260  << "Window Savitzky Golay " << window_savitzky_golay;
261  MOFEM_LOG("WORLD", Sev::inform) << "Window shift " << window_shift;
262
263  if ((rZ.size() - 1) % 2)
264  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
265  "Expected even number of images");
266
267  focalIndex = (rZ.size() - 1) / 2 + window_shift;
268  MOFEM_LOG("WORLD", Sev::inform) << "zR for mid plane " << rZ[focalIndex];
269
270  const auto d1_sg_data = d1_savitzky_golay[order_savitzky_golay];
271  if (!d1_sg_data)
272  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
273  "Wrong Savitzky Golay order");
274
275  auto d1_sg_data_window = d1_sg_data[window_savitzky_golay];
276  if (!d1_sg_data_window)
277  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
278  "Wrong Savitzky Golay window");
279
280  savitzkyGolayNormalisation =
282  savitzkyGolayWeights = d1_sg_data_window;
283
284  CHKERR assembleSystemIntensity();
285  CHKERR solveSystem();
286  CHKERR outputResults(0);
287
289 }
290 //! [run problem]
291
295
296  auto get_bounding_box = [&]() {
297  auto &moab = mField.get_moab();
298  Range verts;
299  MOAB_THROW(moab.get_entities_by_type(0, MBVERTEX, verts));
300
301  ParallelComm *pcomm =
302  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
303
304  Range verts_part;
305  CHKERR pcomm->filter_pstatus(verts, PSTATUS_SHARED | PSTATUS_MULTISHARED,
306  PSTATUS_NOT, -1, &verts_part);
307
308  MatrixDouble coords(verts_part.size(), 3);
309  CHKERR moab.get_coords(verts_part, &*coords.data().begin());
310
311  std::array<double, 2> ave_coords{0, 0};
312  for (auto v = 0; v != verts_part.size(); ++v) {
313  ave_coords[0] += coords(v, 0);
314  ave_coords[1] += coords(v, 1);
315  }
316
317  auto comm = mField.get_comm();
318
319  int local_count = verts_part.size();
320  int global_count;
321  MPI_Allreduce(&local_count, &global_count, 1, MPI_INT, MPI_SUM, comm);
322  std::array<double, 2> ave_coords_glob{0, 0};
323  MPI_Allreduce(ave_coords.data(), ave_coords.data(), 2, MPI_DOUBLE, MPI_SUM,
324  comm);
325  ave_coords_glob[0] /= global_count;
326  ave_coords_glob[1] /= global_count;
327
328  std::array<double, 2> max_coords{ave_coords_glob[0], ave_coords_glob[1]};
329  for (auto v = 0; v != verts_part.size(); ++v) {
330  max_coords[0] = std::max(max_coords[0], coords(v, 0));
331  max_coords[1] = std::max(max_coords[1], coords(v, 1));
332  }
333  std::array<double, 2> max_coords_glob{0, 0};
334  MPI_Allreduce(max_coords.data(), max_coords_glob.data(), 2, MPI_DOUBLE,
335  MPI_MAX, comm);
336
337  std::array<double, 2> min_coords{max_coords_glob[0], max_coords_glob[1]};
338  for (auto v = 0; v != verts_part.size(); ++v) {
339  min_coords[0] = std::min(min_coords[0], coords(v, 0));
340  min_coords[1] = std::min(min_coords[1], coords(v, 1));
341  }
342  std::array<double, 2> min_coords_glob{0, 0};
343  MPI_Allreduce(min_coords.data(), min_coords_glob.data(), 2, MPI_DOUBLE,
344  MPI_MIN, comm);
345
346  return std::array<double, LAST_BB>{ave_coords_glob[0], ave_coords_glob[1],
347  max_coords_glob[0], max_coords_glob[1],
348  min_coords_glob[0], min_coords_glob[1]};
349  };
350
351  CHKERR mField.getInterface(simpleInterface);
352  CHKERR simpleInterface->getOptions();
354
355  aveMaxMin = get_bounding_box();
356
357  MOFEM_LOG("WORLD", Sev::inform)
358  << "Centre " << aveMaxMin[CENTER_X] << " " << aveMaxMin[CENTER_Y];
359  MOFEM_LOG("WORLD", Sev::inform)
360  << "Max " << aveMaxMin[MAX_X] << " " << aveMaxMin[MAX_Y];
361  MOFEM_LOG("WORLD", Sev::inform)
362  << "Min " << aveMaxMin[MIN_X] << " " << aveMaxMin[MIN_Y];
363
365 }
367
368 //! [Set up problem]
371
372  // Note that in 2D case HDIV and HCURL spaces are isomorphic, and therefore
373  // only base for HCURL has been implemented in 2D. Base vectors for HDIV space
374  // are be obtained after rotation of HCURL base vectors by a right angle
375  CHKERR simpleInterface->addDomainField("S", HCURL, DEMKOWICZ_JACOBI_BASE, 1);
377  1);
378  // We use AINSWORTH_LEGENDRE_BASE since DEMKOWICZ_JACOBI_BASE for triangle
379  // is not yet implemented for L2 space. For quads DEMKOWICZ_JACOBI_BASE and
380  // AINSWORTH_LEGENDRE_BASE are construcreed in the same way
381  CHKERR simpleInterface->addDomainField("PHI", L2, DEMKOWICZ_JACOBI_BASE, 1);
382
383  int base_order = 1;
384  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-base_order", &base_order,
385  PETSC_NULL);
386
387
388  MOFEM_LOG("WORLD", Sev::inform) << "Base order " << base_order;
389
390  CHKERR simpleInterface->setFieldOrder("S", base_order);
391  CHKERR simpleInterface->setFieldOrder("PHI", base_order - 1);
392  CHKERR simpleInterface->setUp();
393
394  // Block PHI-PHI is empty
396
398 }
399 //! [Set up problem]
400
401 //! [Calculate flux on boundary]
404  auto pipeline_mng = mField.getInterface<PipelineManager>();
405
406  pipeline_mng->getDomainLhsFE().reset();
407  pipeline_mng->getDomainRhsFE().reset();
408  pipeline_mng->getBoundaryRhsFE().reset();
409
410  auto rule_vol = [](int, int, int order) { return 2 * (order + 1); };
411  pipeline_mng->setBoundaryRhsIntegrationRule(rule_vol);
412
413  auto flux_ptr = boost::make_shared<MatrixDouble>();
414  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
416  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
417  new OpCalculateHVecVectorField<3>("S", flux_ptr));
418  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
419  new BoundaryOp(flux_ptr, calc_flux));
420
421  calc_flux = 0;
422  CHKERR pipeline_mng->loopFiniteElements();
423  double global_flux_assembeld = 0;
424  MPI_Allreduce(&calc_flux, &global_flux_assembeld, 1, MPI_DOUBLE, MPI_SUM,
425  mField.get_comm());
426  calc_flux = global_flux_assembeld;
427
429 }
430 //! [Calculate flux on boundary]
431
432 //! [Push operators to pipeline]
435
436  auto *pipeline_mng = mField.getInterface<PipelineManager>();
437
438  pipeline_mng->getDomainLhsFE().reset();
439  pipeline_mng->getDomainRhsFE().reset();
440  pipeline_mng->getBoundaryRhsFE().reset();
441
442  auto rule_vol = [](int, int, int order) { return 2 * (order + 1); };
443  pipeline_mng->setDomainLhsIntegrationRule(rule_vol);
444  pipeline_mng->setDomainRhsIntegrationRule(rule_vol);
445
446  auto det_ptr = boost::make_shared<VectorDouble>();
447  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
448  auto jac_ptr = boost::make_shared<MatrixDouble>();
449  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
450  pipeline_mng->getOpDomainLhsPipeline().push_back(
451  new OpCalculateHOJacForFace(jac_ptr));
452  pipeline_mng->getOpDomainLhsPipeline().push_back(
453  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
454  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMakeHdivFromHcurl());
455  pipeline_mng->getOpDomainLhsPipeline().push_back(
457  pipeline_mng->getOpDomainLhsPipeline().push_back(
458  new OpSetInvJacHcurlFace(inv_jac_ptr));
459  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
460
461  pipeline_mng->getOpDomainLhsPipeline().push_back(
462  new OpHdivHdiv("S", "S", lhsFlux));
463  auto unity = []() { return 1; };
464  pipeline_mng->getOpDomainLhsPipeline().push_back(
465  new OpHdivU("S", "PHI", unity, true));
466
467  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetHOWeightsOnFace());
468  pipeline_mng->getOpDomainRhsPipeline().push_back(
469  new OpDomainSource("PHI", rhsSource));
470
472 }
473 //! [Push operators to pipeline]
474
475 //! [Solve]
478
479  auto dm = simpleInterface->getDM();
480  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
481
482  auto solver = pipeline_mng->createKSP();
483  CHKERR KSPSetFromOptions(solver);
484  CHKERR KSPSetUp(solver);
485
486  auto iD = createDMVector(dm);
487  auto F = vectorDuplicate(iD);
488  CHKERR KSPSolve(solver, F, iD);
489  CHKERR VecGhostUpdateBegin(iD, INSERT_VALUES, SCATTER_FORWARD);
490  CHKERR VecGhostUpdateEnd(iD, INSERT_VALUES, SCATTER_FORWARD);
491
492  CHKERR DMoFEMMeshToLocalVector(dm, iD, INSERT_VALUES, SCATTER_REVERSE);
493  double i_lambda_flux;
494  CHKERR calculateFlux(i_lambda_flux);
495  MOFEM_LOG_CHANNEL("WORLD");
496  MOFEM_LOG("WORLD", Sev::inform)
497  << "iD flux " << std::scientific << i_lambda_flux;
498
500 }
501 //! [Solve]
502
503 //! [Postprocess results]
506  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
507  pipeline_mng->getDomainLhsFE().reset();
508  pipeline_mng->getDomainRhsFE().reset();
509  pipeline_mng->getBoundaryRhsFE().reset();
510
512
513  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
514  auto jac_ptr = boost::make_shared<MatrixDouble>();
515  post_proc_fe->getOpPtrVector().push_back(
516  new OpCalculateHOJacForFace(jac_ptr));
517  post_proc_fe->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
518  post_proc_fe->getOpPtrVector().push_back(
520
521  auto s_ptr = boost::make_shared<VectorDouble>();
522  auto phi_ptr = boost::make_shared<MatrixDouble>();
523  post_proc_fe->getOpPtrVector().push_back(
524  new OpCalculateScalarFieldValues("S", s_ptr));
525  post_proc_fe->getOpPtrVector().push_back(
526  new OpCalculateHVecVectorField<3>("PHI", phi_ptr));
527
529
530  post_proc_fe->getOpPtrVector().push_back(
531
532  new OpPPMap(post_proc_fe->getPostProcMesh(),
533  post_proc_fe->getMapGaussPts(),
534
535  OpPPMap::DataMapVec{{"S", s_ptr}},
536
537  OpPPMap::DataMapMat{{"PHI", phi_ptr}},
538
540
542
543  )
544
545  );
546
549
550
551  pipeline_mng->getDomainRhsFE() = post_proc_fe;
552  CHKERR pipeline_mng->loopFiniteElements();
553  CHKERR post_proc_fe->writeFile("out_" + boost::lexical_cast<std::string>(i) +
554  ".h5m");
556 }
557 //! [Postprocess results]
558
559 int main(int argc, char *argv[]) {
560
561  // Initialisation of MoFEM/PETSc and MOAB data structures
562  const char param_file[] = "param_file.petsc";
563  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
564
565  try {
566
567  //! [Register MoFEM discrete manager in PETSc]
568  DMType dm_name = "DMMOFEM";
569  CHKERR DMRegister_MoFEM(dm_name);
570  //! [Register MoFEM discrete manager in PETSc
571
572  //! [Create MoAB]
573  moab::Core mb_instance; ///< mesh database
574  moab::Interface &moab = mb_instance; ///< mesh database interface
575  //! [Create MoAB]
576
577  //! [Create MoFEM]
578  MoFEM::Core core(moab); ///< finite element database
579  MoFEM::Interface &m_field = core; ///< finite element database insterface
580  //! [Create MoFEM]
581
582  //! [Example]
583  Example ex(m_field);
584  CHKERR ex.runProblem();
585  //! [Example]
586  }
587  CATCH_ERRORS;
588
590 }
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
MoFEM::EdgeElementForcesAndSourcesCore
Edge finite element.
Definition: EdgeElementForcesAndSourcesCore.hpp:30
MoFEM::EntitiesFieldData::EntData
Data on single entity (This is passed as argument to DataOperator::doWork)
Definition: EntitiesFieldData.hpp:128
Example::lhsFlux
static double lhsFlux(const double x, const double y, const double)
Definition: phase.cpp:165
help
static char help[]
Definition: phase.cpp:17
d1_savitzky_golay_w7_p2
constexpr std::array< int, 7 > d1_savitzky_golay_w7_p2
Definition: phase.cpp:28
Example::MAX_X
@ MAX_X
Definition: phase.cpp:101
Example::aveMaxMin
static std::array< double, LAST_BB > aveMaxMin
Definition: phase.cpp:110
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::OpSetInvJacHcurlFace
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
Definition: UserDataOperators.hpp:2978
Example::CENTER_Y
@ CENTER_Y
Definition: phase.cpp:100
Example::rZ
static std::vector< double > rZ
Definition: phase.cpp:108
OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]
Definition: thermo_elastic.cpp:47
Example::Example
Example(MoFEM::Interface &m_field)
Definition: phase.cpp:81
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
d1_savitzky_golay_w7_p4
constexpr std::array< int, 7 > d1_savitzky_golay_w7_p4
Definition: phase.cpp:33
MoFEM::OpSetContravariantPiolaTransformOnFace2D
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
Definition: UserDataOperators.hpp:3076
Example::BoundaryOp::globFlux
double & globFlux
Definition: phase.cpp:192
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:700
L2
@ L2
field with C-1 continuity
Definition: definitions.h:88
MoFEM::PipelineManager::loopFiniteElements
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Definition: PipelineManager.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
main
int main(int argc, char *argv[])
[Postprocess results]
Definition: phase.cpp:559
d1_savitzky_golay_p4
const int * d1_savitzky_golay_p4[]
Definition: phase.cpp:49
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
d1_savitzky_golay_w9_p2
constexpr std::array< int, 9 > d1_savitzky_golay_w9_p2
Definition: phase.cpp:29
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM.hpp
MoFEM::DMoFEMMeshToLocalVector
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:523
BasicFiniteElements.hpp
MoFEM::PipelineManager::getBoundaryRhsFE
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Definition: PipelineManager.hpp:413
MOAB_THROW
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:554
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
OpBase
OpBaseImpl< PETSC, EdgeEleOp > OpBase
k
static double k
Definition: phase.cpp:75
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
THROW_MESSAGE
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:574
sdf.r
int r
Definition: sdf.py:8
order
constexpr int order
Definition: dg_projection.cpp:18
Example::savitzkyGolayNormalisation
static int savitzkyGolayNormalisation
Definition: phase.cpp:113
Example::calculateFlux
MoFEMErrorCode calculateFlux(double &calc_flux)
[Set up problem]
Definition: phase.cpp:402
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
Example::focalIndex
static int focalIndex
Definition: phase.cpp:112
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
Example::MIN_Y
@ MIN_Y
Definition: phase.cpp:104
d1_normalisation_p2
constexpr std::array< int, 10 > d1_normalisation_p2
Definition: phase.cpp:38
MoFEM::OpBaseImpl
Definition: FormsIntegrators.hpp:178
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:23
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
order_savitzky_golay
static int order_savitzky_golay
Definition: phase.cpp:76
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
[Only used with Hooke equation (linear material model)]
Definition: seepage.cpp:57
d1_savitzky_golay_w5_p4
constexpr std::array< int, 5 > d1_savitzky_golay_w5_p4
Definition: phase.cpp:32
MoFEM::OpCalculateHVecVectorField
Get vector field for H-div approximation.
Definition: UserDataOperators.hpp:2115
Example
[Example]
Definition: plastic.cpp:177
Example::BoundaryOp
Definition: phase.cpp:171
convert.type
type
Definition: convert.py:64
Example::assembleSystemIntensity
MoFEMErrorCode assembleSystemIntensity()
[Calculate flux on boundary]
Definition: phase.cpp:433
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:312
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
Example::solveSystem
MoFEMErrorCode solveSystem()
[Solve]
Definition: dynamic_first_order_con_law.cpp:893
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::PipelineManager::createKSP
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: PipelineManager.cpp:49
Example::rhsSource
static double rhsSource(const double x, const double y, const double)
Definition: phase.cpp:150
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
Example::savitzkyGolayWeights
static const int * savitzkyGolayWeights
Definition: phase.cpp:114
Example::BoundaryOp::doWork
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: phase.cpp:175
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:122
MoFEM::FaceElementForcesAndSourcesCore
Face finite element.
Definition: FaceElementForcesAndSourcesCore.hpp:23
OpHdivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
Definition: phase.cpp:71
Example::getCoordsInImage
static std::pair< int, int > getCoordsInImage(double x, double y)
Definition: phase.cpp:131
Example::BoundaryOp::fluxPtr
boost::shared_ptr< MatrixDouble > fluxPtr
Definition: phase.cpp:191
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
Definition: phase.cpp:73
EntData
EntitiesFieldData::EntData EntData
Definition: phase.cpp:65
d1_normalisation_p4
constexpr std::array< int, 10 > d1_normalisation_p4
Definition: phase.cpp:40
MoFEM::ForcesAndSourcesCore::UserDataOperator
friend class UserDataOperator
Definition: ForcesAndSourcesCore.hpp:482
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:701
BiLinearForm
window_savitzky_golay
static int window_savitzky_golay
Definition: phase.cpp:77
FTensor::Index< 'i', 3 >
MoFEM::EdgeElementForcesAndSourcesCore::UserDataOperator
default operator for EDGE element
Definition: EdgeElementForcesAndSourcesCore.hpp:68
v
const double v
phase velocity of light in medium (cm/ns)
Definition: initial_diffusion.cpp:40
d1_savitzky_golay_p2
const int * d1_savitzky_golay_p2[]
Definition: phase.cpp:43
Range
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
DomainEleOp
MoFEM::CoreTmp< 0 >::Initialize
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
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
OpHdivU
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
Definition: mixed_poisson.cpp:25
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:225
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
Example::BoundaryOp::BoundaryOp
BoundaryOp(boost::shared_ptr< MatrixDouble > flux_ptr, double &glob_flux)
Definition: phase.cpp:172
OpHdivHdiv
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Definition: phase.cpp:69
HCURL
@ HCURL
field with continuous tangents
Definition: definitions.h:86
d1_savitzky_golay
const int ** d1_savitzky_golay[]
Definition: phase.cpp:60
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::PetscOptionsGetString
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
Definition: DeprecatedPetsc.hpp:172
Example::MAX_Y
@ MAX_Y
Definition: phase.cpp:102
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
MoFEM::OpCalculateHOJacForFace
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
Definition: HODataOperators.hpp:264
Example::iI
static std::vector< MatrixInt > iI
Definition: phase.cpp:109
Example::runProblem
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:213
d1_savitzky_golay_w5_p2
constexpr std::array< int, 5 > d1_savitzky_golay_w5_p2
Definition: phase.cpp:27
sdf_wavy_2d.w
int w
Definition: sdf_wavy_2d.py:6
m
FTensor::Index< 'm', 3 > m
Definition: shallow_wave.cpp:80
d1_normalisation
const int * d1_normalisation[]
Definition: phase.cpp:62
MoFEM::PetscOptionsGetScalar
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:162
MoFEM::OpMakeHdivFromHcurl
Make Hdiv space from Hcurl space in 2d.
Definition: UserDataOperators.hpp:2985
DomainEle
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Definition: child_and_parent.cpp:34
OpDomainSource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
Definition: child_and_parent.cpp:55
Example::MIN_X
@ MIN_X
Definition: phase.cpp:103
Example::outputResults
MoFEMErrorCode outputResults()
[Solve]
Definition: dynamic_first_order_con_law.cpp:1183
convert.int
int
Definition: convert.py:64
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
d1_savitzky_golay_w9_p4
constexpr std::array< int, 9 > d1_savitzky_golay_w9_p4
Definition: phase.cpp:35
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
Example::BoundingBox
BoundingBox
Definition: phase.cpp:98
d1_savitzky_golay_w3_p2
constexpr std::array< int, 3 > d1_savitzky_golay_w3_p2
Definition: phase.cpp:26
F
@ F
Definition: free_surface.cpp:394
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698