v0.13.1
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
15using namespace MoFEM;
16
17static char help[] = "...\n\n";
18
19#include <BasicFiniteElements.hpp>
20
25
26constexpr std::array<int, 3> d1_savitzky_golay_w3_p2 = {-1, 0, 1};
27constexpr std::array<int, 5> d1_savitzky_golay_w5_p2 = {-2, -1, 0, 1, 2};
28constexpr std::array<int, 7> d1_savitzky_golay_w7_p2 = {-3, -2, -1, 0, 1, 2, 3};
29constexpr std::array<int, 9> d1_savitzky_golay_w9_p2 = {-4, -3, -2, -1, 0,
30 1, 2, 3, 4};
31
32constexpr std::array<int, 5> d1_savitzky_golay_w5_p4 = {1, -8, 0, 8, -1};
33constexpr std::array<int, 7> d1_savitzky_golay_w7_p4 = {22, -67, -58, 0,
34 58, 67, -22};
35constexpr std::array<int, 9> d1_savitzky_golay_w9_p4 = {
36 86, -142, -193, -126, 0, 126, 193, 142, -86};
37
38constexpr std::array<int, 10> d1_normalisation_p2 = {0, 0, 0, 2, 0,
39 10, 0, 28, 0, 60};
40constexpr std::array<int, 10> d1_normalisation_p4 = {0, 0, 0, 0, 0,
41 12, 0, 252, 0, 1188};
42
43const 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
49const int *d1_savitzky_golay_p4[] = {nullptr,
50 nullptr,
51 nullptr,
52 nullptr,
53 nullptr,
55 nullptr,
57 nullptr,
59
60const int **d1_savitzky_golay[] = {nullptr, nullptr, d1_savitzky_golay_p2,
61 nullptr, d1_savitzky_golay_p4};
62const int *d1_normalisation[] = {nullptr, nullptr, d1_normalisation_p2.data(),
63 nullptr, d1_normalisation_p4.data()};
64
71 GAUSS>::OpMixDivTimesScalar<2>;
74
75static double k = 1; // wave number
76static int order_savitzky_golay = 2;
77static int window_savitzky_golay = 3;
78
79struct Example {
80
81 Example(MoFEM::Interface &m_field) : mField(m_field) {}
82
84
85private:
88 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
89
95 MoFEMErrorCode calculateFlux(double &calc_flux);
97
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
124std::vector<double> Example::rZ;
125std::vector<MatrixInt> Example::iI;
126std::array<double, Example::BoundingBox::LAST_BB> Example::aveMaxMin;
130
131std::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
150double 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
165double 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
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
190private:
191 boost::shared_ptr<MatrixDouble> fluxPtr;
192 double &globFlux;
193};
194
195//! [run problem]
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
282 savitzkyGolayWeights = d1_sg_data_window;
283
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
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
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
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);
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();
481
482 auto solver = pipeline_mng->createKSP();
483 CHKERR KSPSetFromOptions(solver);
484 CHKERR KSPSetUp(solver);
485
486 auto iD = smartCreateDMVector(dm);
487 auto F = smartVectorDuplicate(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]
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
559int main(int argc, char *argv[]) {
560
561 // Initialisation of MoFEM/PETSc and MOAB data structures
562 const char param_file[] = "param_file.petsc";
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 }
588
590}
static Index< 'p', 3 > p
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
Definition: definitions.h:541
@ L2
field with C-1 continuity
Definition: definitions.h:88
@ NOSPACE
Definition: definitions.h:83
@ HCURL
field with continuous tangents
Definition: definitions.h:86
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:215
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#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
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
Definition: definitions.h:561
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
FTensor::Index< 'm', SPACE_DIM > m
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:301
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:277
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
OpSetInvJacHcurlFaceImpl< 2 > OpSetInvJacHcurlFace
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
OpSetContravariantPiolaTransformOnFace2DImpl< 2 > OpSetContravariantPiolaTransformOnFace2D
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
double w(const double g, const double t)
const double r
rate factor
int main(int argc, char *argv[])
[Postprocess results]
Definition: phase.cpp:559
const int * d1_normalisation[]
Definition: phase.cpp:62
constexpr std::array< int, 3 > d1_savitzky_golay_w3_p2
Definition: phase.cpp:26
EntitiesFieldData::EntData EntData
Definition: phase.cpp:65
constexpr std::array< int, 10 > d1_normalisation_p2
Definition: phase.cpp:38
static char help[]
Definition: phase.cpp:17
constexpr std::array< int, 9 > d1_savitzky_golay_w9_p4
Definition: phase.cpp:35
constexpr std::array< int, 7 > d1_savitzky_golay_w7_p4
Definition: phase.cpp:33
constexpr std::array< int, 5 > d1_savitzky_golay_w5_p2
Definition: phase.cpp:27
constexpr std::array< int, 10 > d1_normalisation_p4
Definition: phase.cpp:40
static int window_savitzky_golay
Definition: phase.cpp:77
static int order_savitzky_golay
Definition: phase.cpp:76
const int ** d1_savitzky_golay[]
Definition: phase.cpp:60
constexpr std::array< int, 7 > d1_savitzky_golay_w7_p2
Definition: phase.cpp:28
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
Definition: phase.cpp:71
constexpr std::array< int, 5 > d1_savitzky_golay_w5_p4
Definition: phase.cpp:32
const int * d1_savitzky_golay_p2[]
Definition: phase.cpp:43
const int * d1_savitzky_golay_p4[]
Definition: phase.cpp:49
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Definition: phase.cpp:69
static double k
Definition: phase.cpp:75
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpDomainSource
Definition: phase.cpp:73
constexpr std::array< int, 9 > d1_savitzky_golay_w9_p2
Definition: phase.cpp:29
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
double scale
Definition: plastic.cpp:85
OpBaseImpl< PETSC, EdgeEleOp > OpBase
boost::shared_ptr< MatrixDouble > fluxPtr
Definition: phase.cpp:191
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Definition: phase.cpp:175
double & globFlux
Definition: phase.cpp:192
BoundaryOp(boost::shared_ptr< MatrixDouble > flux_ptr, double &glob_flux)
Definition: phase.cpp:172
[Example]
Definition: plastic.cpp:111
static std::array< double, LAST_BB > aveMaxMin
Definition: phase.cpp:110
static int focalIndex
Definition: phase.cpp:112
static std::vector< MatrixInt > iI
Definition: phase.cpp:109
static double rhsSource(const double x, const double y, const double)
Definition: phase.cpp:150
MoFEMErrorCode solveSystem()
static std::pair< int, int > getCoordsInImage(double x, double y)
Definition: phase.cpp:131
Example(MoFEM::Interface &m_field)
Definition: phase.cpp:81
Simple * simpleInterface
Definition: helmholtz.cpp:45
MoFEMErrorCode calculateFlux(double &calc_flux)
[Set up problem]
Definition: phase.cpp:402
MoFEMErrorCode runProblem()
MoFEMErrorCode assembleSystemFlux()
MoFEM::Interface & mField
Definition: plastic.cpp:118
MoFEMErrorCode assembleSystemIntensity()
[Calculate flux on boundary]
Definition: phase.cpp:433
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:137
static double lhsFlux(const double x, const double y, const double)
Definition: phase.cpp:165
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
[Solve]
Definition: helmholtz.cpp:255
BoundingBox
Definition: phase.cpp:98
@ MAX_X
Definition: phase.cpp:101
@ MIN_X
Definition: phase.cpp:103
@ MIN_Y
Definition: phase.cpp:104
@ CENTER_X
Definition: phase.cpp:99
@ MAX_Y
Definition: phase.cpp:102
@ CENTER_Y
Definition: phase.cpp:100
@ LAST_BB
Definition: phase.cpp:105
static int savitzkyGolayNormalisation
Definition: phase.cpp:113
static std::vector< double > rZ
Definition: phase.cpp:108
static const int * savitzkyGolayWeights
Definition: phase.cpp:114
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =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.
auto getFTensor1Normal()
get edge normal NOTE: it should be used only in 2D analysis
Data on single entity (This is passed as argument to DataOperator::doWork)
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPSPACE
operator do Work is execute on space data
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get vector field for H-div approximation.
Get value at integration points for scalar field.
Make Hdiv space from Hcurl space in 2d.
Post post-proc data at points from hash maps.
Definition: PostProc.hpp:166
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProc.hpp:168
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProc.hpp:169
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Definition: PostProc.hpp:120
Simple interface for fast problem set-up.
Definition: Simple.hpp:26
MoFEMErrorCode addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition: Simple.cpp:374
MoFEMErrorCode addFieldToEmptyFieldBlocks(const std::string row_field, const std::string col_field) const
add empty block to problem
Definition: Simple.cpp:841
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode addBoundaryField(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 boundary.
Definition: Simple.cpp:392
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, 3 > OpHdivHdiv
Integrate Lhs base of flux (1/k) base of flux (FLUX x FLUX)