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 
21 using DomainEle = FaceElementForcesAndSourcesCore;
22 using DomainEleOp = DomainEle::UserDataOperator;
23 using EdgeEle = EdgeElementForcesAndSourcesCore;
24 using EdgeEleOp = EdgeEle::UserDataOperator;
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,
54  d1_savitzky_golay_w5_p4.data(),
55  nullptr,
56  d1_savitzky_golay_w7_p4.data(),
57  nullptr,
58  d1_savitzky_golay_w9_p4.data()};
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 
65 using EntData = EntitiesFieldData::EntData;
66 using AssemblyBoundaryEleOp =
67  FormsIntegrators<EdgeEleOp>::Assembly<PETSC>::OpBase;
68 using OpHdivHdiv = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
69  GAUSS>::OpMass<3, 3>;
70 using OpHdivU = FormsIntegrators<DomainEleOp>::Assembly<PETSC>::BiLinearForm<
71  GAUSS>::OpMixDivTimesScalar<2>;
72 using OpDomainSource = FormsIntegrators<DomainEleOp>::Assembly<
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 
90  MoFEMErrorCode readMesh();
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,
100  CENTER_Y,
101  MAX_X,
102  MAX_Y,
103  MIN_X,
104  MIN_Y,
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;
113  static int savitzkyGolayNormalisation;
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;
127 int Example::focalIndex;
128 int Example::savitzkyGolayNormalisation;
129 const int *Example::savitzkyGolayWeights;
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 
171 struct Example::BoundaryOp : public EdgeEleOp {
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) {
176  MoFEMFunctionBegin;
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();
181  FTensor::Index<'i', 3> i;
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  };
187  MoFEMFunctionReturn(0);
188  }
189 
190 private:
191  boost::shared_ptr<MatrixDouble> fluxPtr;
192  double &globFlux;
193 };
194 
195 //! [run problem]
196 MoFEMErrorCode Example::runProblem() {
197  MoFEMFunctionBegin;
198  CHKERR readMesh();
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 
244  structure_data(read_images());
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 =
281  d1_normalisation[order_savitzky_golay][window_savitzky_golay];
282  savitzkyGolayWeights = d1_sg_data_window;
283 
284  CHKERR assembleSystemIntensity();
285  CHKERR solveSystem();
286  CHKERR outputResults(0);
287 
288  MoFEMFunctionReturn(0);
289 }
290 //! [run problem]
291 
292 //! [Read mesh]
293 MoFEMErrorCode Example::readMesh() {
294  MoFEMFunctionBegin;
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();
353  CHKERR simpleInterface->loadFile();
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 
364  MoFEMFunctionReturn(0);
365 }
366 //! [Read mesh]
367 
368 //! [Set up problem]
369 MoFEMErrorCode Example::setupProblem() {
370  MoFEMFunctionBegin;
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);
376  CHKERR simpleInterface->addBoundaryField("S", HCURL, DEMKOWICZ_JACOBI_BASE,
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
395  CHKERR simpleInterface->addFieldToEmptyFieldBlocks("PHI", "PHI");
396 
397  MoFEMFunctionReturn(0);
398 }
399 //! [Set up problem]
400 
401 //! [Calculate flux on boundary]
402 MoFEMErrorCode Example::calculateFlux(double &calc_flux) {
403  MoFEMFunctionBegin;
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(
415  new OpSetContravariantPiolaTransformOnEdge2D());
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 
428  MoFEMFunctionReturn(0);
429 }
430 //! [Calculate flux on boundary]
431 
432 //! [Push operators to pipeline]
433 MoFEMErrorCode Example::assembleSystemIntensity() {
434  MoFEMFunctionBegin;
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(
456  new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
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 
471  MoFEMFunctionReturn(0);
472 }
473 //! [Push operators to pipeline]
474 
475 //! [Solve]
476 MoFEMErrorCode Example::solveSystem() {
477  MoFEMFunctionBegin;
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 
499  MoFEMFunctionReturn(0);
500 }
501 //! [Solve]
502 
503 //! [Postprocess results]
504 MoFEMErrorCode Example::outputResults(const int i) {
505  MoFEMFunctionBegin;
506  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
507  pipeline_mng->getDomainLhsFE().reset();
508  pipeline_mng->getDomainRhsFE().reset();
509  pipeline_mng->getBoundaryRhsFE().reset();
510 
511  using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
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(
519  new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
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 
528  using OpPPMap = OpPostProcMapInMoab<3, 3>;
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 
539  OpPPMap::DataMapMat{},
540 
541  OpPPMap::DataMapMat{}
542 
543  )
544 
545  );
546 
547  // post_proc_fe->addFieldValuesPostProc("S");
548  // post_proc_fe->addFieldValuesPostProc("PHI");
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");
555  MoFEMFunctionReturn(0);
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 
589  CHKERR MoFEM::Core::Finalize();
590 }
NOSPACE
@ NOSPACE
Definition: definitions.h:83
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce 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:127
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:215
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
MoFEM::OpSetInvJacHcurlFace
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
Definition: UserDataOperators.hpp:2982
PlasticOps::w
double w(double eqiv, double dot_tau, double f, double sigma_y, double sigma_Y)
Definition: PlasticOpsGeneric.hpp:276
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:3080
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:104
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:527
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:541
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
Definition: radiation.cpp:29
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:561
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:1975
Example::focalIndex
static int focalIndex
Definition: phase.cpp:112
Example::readMesh
MoFEMErrorCode readMesh()
[Run problem]
Definition: dynamic_first_order_con_law.cpp:463
c
const double c
speed of light (cm/ns)
Definition: initial_diffusion.cpp:39
MoFEM::PostProcBrokenMeshInMoab
Definition: PostProcBrokenMeshInMoabBase.hpp:667
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:170
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1018
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:2116
Example
[Example]
Definition: plastic.cpp:228
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:302
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::OpSetContravariantPiolaTransformOnEdge2D
Definition: UserDataOperators.hpp:3090
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:47
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:128
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:217
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
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:1094
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:276
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:3254
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:264
d1_savitzky_golay_w5_p2
constexpr std::array< int, 5 > d1_savitzky_golay_w5_p2
Definition: phase.cpp:27
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:2989
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:416
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:346
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
Generated by Doxygen 1.8.17 and hosted at University of Glasgow