v0.14.0
phase.cpp

Calculate light phase using Transport intensity equation.

See paper [55]

/**
* \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";
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>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
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) {}
MoFEMErrorCode runProblem();
private:
Simple *simpleInterface;
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
MoFEMErrorCode readMesh();
MoFEMErrorCode setupProblem();
MoFEMErrorCode assembleSystemIntensity();
MoFEMErrorCode assembleSystemFlux();
MoFEMErrorCode solveSystem();
MoFEMErrorCode calculateFlux(double &calc_flux);
MoFEMErrorCode outputResults(const int i);
enum BoundingBox {
CENTER_X = 0,
CENTER_Y,
MAX_X,
MAX_Y,
MIN_X,
MIN_Y,
LAST_BB
};
static std::vector<double> rZ;
static std::vector<MatrixInt> iI;
static std::array<double, LAST_BB> aveMaxMin;
static int focalIndex;
static int savitzkyGolayNormalisation;
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 -= aveMaxMin[MIN_X];
y -= aveMaxMin[MIN_Y];
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) {}
MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
const auto nb_gauss_pts = getGaussPts().size2();
auto t_flux = getFTensor1FromMat<3>(*fluxPtr);
auto t_normal = getFTensor1Normal();
auto t_w = getFTensor0IntegrationWeight();
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]
CHKERR readMesh();
CHKERR setupProblem();
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++);
}
}
}
};
structure_data(read_images());
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");
savitzkyGolayNormalisation =
savitzkyGolayWeights = d1_sg_data_window;
CHKERR assembleSystemIntensity();
CHKERR solveSystem();
CHKERR outputResults(0);
}
//! [run problem]
//! [Read mesh]
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 mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
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];
}
//! [Read mesh]
//! [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
CHKERR simpleInterface->addDomainField("S", HCURL, DEMKOWICZ_JACOBI_BASE, 1);
CHKERR simpleInterface->addBoundaryField("S", HCURL, DEMKOWICZ_JACOBI_BASE,
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);
CHKERR simpleInterface->setUp();
// Block PHI-PHI is empty
CHKERR simpleInterface->addFieldToEmptyFieldBlocks("PHI", "PHI");
}
//! [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,
mField.get_comm());
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();
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
CHKERR KSPSetUp(solver);
auto iD = createDMVector(dm);
auto F = vectorDuplicate(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]
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
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}},
)
);
// post_proc_fe->addFieldValuesPostProc("S");
// post_proc_fe->addFieldValuesPostProc("PHI");
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";
MoFEM::Core::Initialize(&argc, &argv, param_file, help);
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]
}
}
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: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::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
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
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
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: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
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:535
ContactOps::scale
double scale
Definition: EshelbianContact.hpp:22
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1067
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: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:310
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:43
Example::savitzkyGolayWeights
static const int * savitzkyGolayWeights
Definition: phase.cpp:114
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
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: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:1102
Example::setupProblem
MoFEMErrorCode setupProblem()
[Run problem]
Definition: plastic.cpp:225
UserDataOperator
ForcesAndSourcesCore::UserDataOperator UserDataOperator
Definition: HookeElement.hpp:75
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
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: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:7
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::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
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