v0.13.1
phase.cpp

Calculate light phase using Transport intensity equation.

See paper [38]

/**
* \file phase.cpp
* \example phase.cpp
*
* Calculate light phase using <a
* href="https://en.wikipedia.org/wiki/Transport-of-intensity_equation">Transport
* intensity equation</a>.
*
* See paper \cite parvizi2015practical
*
*/
#include <MoFEM.hpp>
using namespace MoFEM;
static char help[] = "...\n\n";
#include <BasicFiniteElements.hpp>
constexpr std::array<int, 3> d1_savitzky_golay_w3_p2 = {-1, 0, 1};
constexpr std::array<int, 5> d1_savitzky_golay_w5_p2 = {-2, -1, 0, 1, 2};
constexpr std::array<int, 7> d1_savitzky_golay_w7_p2 = {-3, -2, -1, 0, 1, 2, 3};
constexpr std::array<int, 9> d1_savitzky_golay_w9_p2 = {-4, -3, -2, -1, 0,
1, 2, 3, 4};
constexpr std::array<int, 5> d1_savitzky_golay_w5_p4 = {1, -8, 0, 8, -1};
constexpr std::array<int, 7> d1_savitzky_golay_w7_p4 = {22, -67, -58, 0,
58, 67, -22};
constexpr std::array<int, 9> d1_savitzky_golay_w9_p4 = {
86, -142, -193, -126, 0, 126, 193, 142, -86};
constexpr std::array<int, 10> d1_normalisation_p2 = {0, 0, 0, 2, 0,
10, 0, 28, 0, 60};
constexpr std::array<int, 10> d1_normalisation_p4 = {0, 0, 0, 0, 0,
12, 0, 252, 0, 1188};
const int *d1_savitzky_golay_p2[] = {nullptr, nullptr,
nullptr, d1_savitzky_golay_w3_p2.data(),
nullptr, d1_savitzky_golay_w5_p2.data(),
nullptr, d1_savitzky_golay_w7_p2.data(),
nullptr, d1_savitzky_golay_w9_p2.data()};
const int *d1_savitzky_golay_p4[] = {nullptr,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr,
const int **d1_savitzky_golay[] = {nullptr, nullptr, d1_savitzky_golay_p2,
const int *d1_normalisation[] = {nullptr, nullptr, d1_normalisation_p2.data(),
nullptr, d1_normalisation_p4.data()};
GAUSS>::OpMixDivTimesScalar<2>;
static double k = 1; // wave number
static int order_savitzky_golay = 2;
static int window_savitzky_golay = 3;
struct Example {
Example(MoFEM::Interface &m_field) : mField(m_field) {}
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
MoFEMErrorCode calculateFlux(double &calc_flux);
enum BoundingBox {
CENTER_X = 0,
};
static std::vector<double> rZ;
static std::vector<MatrixInt> iI;
static std::array<double, LAST_BB> aveMaxMin;
static int focalIndex;
static const int *savitzkyGolayWeights;
static std::pair<int, int> getCoordsInImage(double x, double y);
static double rhsSource(const double x, const double y, const double);
static double lhsFlux(const double x, const double y, const double);
struct BoundaryOp;
};
std::vector<double> Example::rZ;
std::vector<MatrixInt> Example::iI;
std::array<double, Example::BoundingBox::LAST_BB> Example::aveMaxMin;
std::pair<int, int> Example::getCoordsInImage(double x, double y) {
auto &m = iI[focalIndex];
x *= (m.size1() - 1) / (aveMaxMin[MAX_X] - aveMaxMin[MIN_X]);
y *= (m.size2() - 1) / (aveMaxMin[MAX_Y] - aveMaxMin[MIN_Y]);
const auto p = std::make_pair<int, int>(std::round(x), std::round(y));
#ifndef NDEBUG
if (p.first < 0 && p.first >= m.size1())
THROW_MESSAGE("Wrong index");
if (p.second < 0 && p.second >= m.size2())
THROW_MESSAGE("Wrong index");
#endif
return p;
}
double Example::rhsSource(const double x, const double y, const double) {
const auto idx = getCoordsInImage(x, y);
double v = 0;
for (auto w = 0; w != window_savitzky_golay; ++w) {
const auto i = focalIndex - (window_savitzky_golay - 1) / 2 + w;
const auto &intensity = iI[i];
v += intensity(idx.first, idx.second) * savitzkyGolayWeights[w];
}
v = static_cast<double>(v) / savitzkyGolayNormalisation;
const auto dz = rZ[focalIndex + 1] - rZ[focalIndex - 1];
return -k * v / dz;
}
double Example::lhsFlux(const double x, const double y, const double) {
const auto idx = getCoordsInImage(x, y);
const auto &m = iI[focalIndex];
return 1. / m(idx.first, idx.second);
}
struct Example::BoundaryOp : public EdgeEleOp {
BoundaryOp(boost::shared_ptr<MatrixDouble> flux_ptr, double &glob_flux)
: EdgeEleOp(NOSPACE, OPSPACE), fluxPtr(flux_ptr), globFlux(glob_flux) {}
const auto nb_gauss_pts = getGaussPts().size2();
auto t_flux = getFTensor1FromMat<3>(*fluxPtr);
auto t_normal = getFTensor1Normal();
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
globFlux += t_w * t_normal(i) * t_flux(i);
++t_flux;
++t_w;
};
}
private:
boost::shared_ptr<MatrixDouble> fluxPtr;
double &globFlux;
};
//! [run problem]
char images_array_files[255] = "out_arrays.txt";
CHKERR PetscOptionsGetString(PETSC_NULL, "", "-images_arrays",
images_array_files, 255, PETSC_NULL);
// Look into "extract_data.ipynb" file, it is used to generate data file with
// combined stack of the images.
auto read_images = [&]() {
std::ifstream in;
in.open(images_array_files);
std::vector<int> values;
values.insert(values.end(), std::istream_iterator<int>(in),
std::istream_iterator<int>());
MOFEM_LOG("WORLD", Sev::inform) << "Read data size " << values.size();
in.close();
return values;
};
auto structure_data = [&](auto &&data) {
constexpr double scale = 1e4; // scale to float
auto it = data.begin();
if (it == data.end()) {
THROW_MESSAGE("No images");
}
rZ.reserve(*it);
iI.reserve(*it);
MOFEM_LOG("WORLD", Sev::inform) << "Number of images " << *it;
it++;
for (; it != data.end();) {
rZ.emplace_back(*(it++) / scale);
const auto r = *(it++);
const auto c = *(it++);
iI.emplace_back(r, c);
MOFEM_LOG("WORLD", Sev::inform)
<< "Read data set " << rZ.back() << " size " << r << " by " << c;
auto &m = iI.back();
for (auto rit = m.begin1(); rit != m.end1(); ++rit) {
for (auto cit = rit.begin(); cit != rit.end(); ++cit) {
*cit = *(it++);
}
}
}
};
CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-k", &k, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order_savitzky_golay",
&order_savitzky_golay, PETSC_NULL);
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-window_savitzky_golay",
&window_savitzky_golay, PETSC_NULL);
int window_shift = 0;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-window_shift", &window_shift,
PETSC_NULL);
MOFEM_LOG("WORLD", Sev::inform) << "Wave number " << k;
MOFEM_LOG("WORLD", Sev::inform)
<< "Order Savitzky Golay " << order_savitzky_golay;
MOFEM_LOG("WORLD", Sev::inform)
<< "Window Savitzky Golay " << window_savitzky_golay;
MOFEM_LOG("WORLD", Sev::inform) << "Window shift " << window_shift;
if ((rZ.size() - 1) % 2)
SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
"Expected even number of images");
focalIndex = (rZ.size() - 1) / 2 + window_shift;
MOFEM_LOG("WORLD", Sev::inform) << "zR for mid plane " << rZ[focalIndex];
const auto d1_sg_data = d1_savitzky_golay[order_savitzky_golay];
if (!d1_sg_data)
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Wrong Savitzky Golay order");
auto d1_sg_data_window = d1_sg_data[window_savitzky_golay];
if (!d1_sg_data_window)
SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
"Wrong Savitzky Golay window");
savitzkyGolayWeights = d1_sg_data_window;
}
//! [run problem]
auto get_bounding_box = [&]() {
auto &moab = mField.get_moab();
Range verts;
MOAB_THROW(moab.get_entities_by_type(0, MBVERTEX, verts));
ParallelComm *pcomm =
ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
Range verts_part;
CHKERR pcomm->filter_pstatus(verts, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &verts_part);
MatrixDouble coords(verts_part.size(), 3);
CHKERR moab.get_coords(verts_part, &*coords.data().begin());
std::array<double, 2> ave_coords{0, 0};
for (auto v = 0; v != verts_part.size(); ++v) {
ave_coords[0] += coords(v, 0);
ave_coords[1] += coords(v, 1);
}
auto comm = mField.get_comm();
int local_count = verts_part.size();
int global_count;
MPI_Allreduce(&local_count, &global_count, 1, MPI_INT, MPI_SUM, comm);
std::array<double, 2> ave_coords_glob{0, 0};
MPI_Allreduce(ave_coords.data(), ave_coords.data(), 2, MPI_DOUBLE, MPI_SUM,
comm);
ave_coords_glob[0] /= global_count;
ave_coords_glob[1] /= global_count;
std::array<double, 2> max_coords{ave_coords_glob[0], ave_coords_glob[1]};
for (auto v = 0; v != verts_part.size(); ++v) {
max_coords[0] = std::max(max_coords[0], coords(v, 0));
max_coords[1] = std::max(max_coords[1], coords(v, 1));
}
std::array<double, 2> max_coords_glob{0, 0};
MPI_Allreduce(max_coords.data(), max_coords_glob.data(), 2, MPI_DOUBLE,
MPI_MAX, comm);
std::array<double, 2> min_coords{max_coords_glob[0], max_coords_glob[1]};
for (auto v = 0; v != verts_part.size(); ++v) {
min_coords[0] = std::min(min_coords[0], coords(v, 0));
min_coords[1] = std::min(min_coords[1], coords(v, 1));
}
std::array<double, 2> min_coords_glob{0, 0};
MPI_Allreduce(min_coords.data(), min_coords_glob.data(), 2, MPI_DOUBLE,
MPI_MIN, comm);
return std::array<double, LAST_BB>{ave_coords_glob[0], ave_coords_glob[1],
max_coords_glob[0], max_coords_glob[1],
min_coords_glob[0], min_coords_glob[1]};
};
CHKERR simpleInterface->getOptions();
aveMaxMin = get_bounding_box();
MOFEM_LOG("WORLD", Sev::inform)
<< "Centre " << aveMaxMin[CENTER_X] << " " << aveMaxMin[CENTER_Y];
MOFEM_LOG("WORLD", Sev::inform)
<< "Max " << aveMaxMin[MAX_X] << " " << aveMaxMin[MAX_Y];
MOFEM_LOG("WORLD", Sev::inform)
<< "Min " << aveMaxMin[MIN_X] << " " << aveMaxMin[MIN_Y];
}
//! [Set up problem]
// Note that in 2D case HDIV and HCURL spaces are isomorphic, and therefore
// only base for HCURL has been implemented in 2D. Base vectors for HDIV space
// are be obtained after rotation of HCURL base vectors by a right angle
1);
// We use AINSWORTH_LEGENDRE_BASE since DEMKOWICZ_JACOBI_BASE for triangle
// is not yet implemented for L2 space. For quads DEMKOWICZ_JACOBI_BASE and
// AINSWORTH_LEGENDRE_BASE are construcreed in the same way
CHKERR simpleInterface->addDomainField("PHI", L2, DEMKOWICZ_JACOBI_BASE, 1);
int base_order = 1;
CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-base_order", &base_order,
PETSC_NULL);
MOFEM_LOG("WORLD", Sev::inform) << "Base order " << base_order;
CHKERR simpleInterface->setFieldOrder("S", base_order);
CHKERR simpleInterface->setFieldOrder("PHI", base_order - 1);
// Block PHI-PHI is empty
}
//! [Set up problem]
//! [Calculate flux on boundary]
auto pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
auto rule_vol = [](int, int, int order) { return 2 * (order + 1); };
pipeline_mng->setBoundaryRhsIntegrationRule(rule_vol);
auto flux_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new OpCalculateHVecVectorField<3>("S", flux_ptr));
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
new BoundaryOp(flux_ptr, calc_flux));
calc_flux = 0;
CHKERR pipeline_mng->loopFiniteElements();
double global_flux_assembeld = 0;
MPI_Allreduce(&calc_flux, &global_flux_assembeld, 1, MPI_DOUBLE, MPI_SUM,
calc_flux = global_flux_assembeld;
}
//! [Calculate flux on boundary]
//! [Push operators to pipeline]
auto *pipeline_mng = mField.getInterface<PipelineManager>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
auto rule_vol = [](int, int, int order) { return 2 * (order + 1); };
pipeline_mng->setDomainLhsIntegrationRule(rule_vol);
pipeline_mng->setDomainRhsIntegrationRule(rule_vol);
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpCalculateHOJacForFace(jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMakeHdivFromHcurl());
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpSetInvJacHcurlFace(inv_jac_ptr));
pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpHdivHdiv("S", "S", lhsFlux));
auto unity = []() { return 1; };
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpHdivU("S", "PHI", unity, true));
pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetHOWeightsOnFace());
pipeline_mng->getOpDomainRhsPipeline().push_back(
new OpDomainSource("PHI", rhsSource));
}
//! [Push operators to pipeline]
//! [Solve]
auto dm = simpleInterface->getDM();
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto iD = smartCreateDMVector(dm);
auto F = smartVectorDuplicate(iD);
CHKERR KSPSolve(solver, F, iD);
CHKERR VecGhostUpdateBegin(iD, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(iD, INSERT_VALUES, SCATTER_FORWARD);
CHKERR DMoFEMMeshToLocalVector(dm, iD, INSERT_VALUES, SCATTER_REVERSE);
double i_lambda_flux;
CHKERR calculateFlux(i_lambda_flux);
MOFEM_LOG("WORLD", Sev::inform)
<< "iD flux " << std::scientific << i_lambda_flux;
}
//! [Solve]
//! [Postprocess results]
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
auto jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHOJacForFace(jac_ptr));
post_proc_fe->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
post_proc_fe->getOpPtrVector().push_back(
auto s_ptr = boost::make_shared<VectorDouble>();
auto phi_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateScalarFieldValues("S", s_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpCalculateHVecVectorField<3>("PHI", phi_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
OpPPMap::DataMapVec{{"S", s_ptr}},
OpPPMap::DataMapMat{{"PHI", phi_ptr}},
OpPPMap::DataMapMat{},
OpPPMap::DataMapMat{}
)
);
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile("out_" + boost::lexical_cast<std::string>(i) +
".h5m");
}
//! [Postprocess results]
int main(int argc, char *argv[]) {
// Initialisation of MoFEM/PETSc and MOAB data structures
const char param_file[] = "param_file.petsc";
try {
//! [Register MoFEM discrete manager in PETSc]
DMType dm_name = "DMMOFEM";
//! [Register MoFEM discrete manager in PETSc
//! [Create MoAB]
moab::Core mb_instance; ///< mesh database
moab::Interface &moab = mb_instance; ///< mesh database interface
//! [Create MoAB]
//! [Create MoFEM]
MoFEM::Core core(moab); ///< finite element database
MoFEM::Interface &m_field = core; ///< finite element database insterface
//! [Create MoFEM]
//! [Example]
Example ex(m_field);
CHKERR ex.runProblem();
//! [Example]
}
}
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
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
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
[run problem]
Definition: helmholtz.cpp:73
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()
[Push operators to pipeline]
Definition: helmholtz.cpp:235
static std::pair< int, int > getCoordsInImage(double x, double y)
Definition: phase.cpp:131
Example(MoFEM::Interface &m_field)
Definition: plastic.cpp:113
Simple * simpleInterface
Definition: helmholtz.cpp:45
MoFEMErrorCode calculateFlux(double &calc_flux)
[Set up problem]
Definition: phase.cpp:402
MoFEMErrorCode runProblem()
[Run problem]
Definition: plastic.cpp:144
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()
[Run problem]
Definition: plastic.cpp:156
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
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 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)