static char help[] =
"...\n\n";
1, 2, 3, 4};
58, 67, -22};
86, -142, -193, -126, 0, 126, 193, 142, -86};
10, 0, 28, 0, 60};
12, 0, 252, 0, 1188};
nullptr,
nullptr,
nullptr,
nullptr,
nullptr,
nullptr,
GAUSS>::OpMixDivTimesScalar<2>;
private:
};
static std::vector<double>
rZ;
static std::vector<MatrixInt>
iI;
static std::array<double, LAST_BB>
aveMaxMin;
static double rhsSource(
const double x,
const double y,
const double);
static double lhsFlux(
const double x,
const double y,
const double);
struct BoundaryOp;
};
const auto p = std::make_pair<int, int>(std::round(x), std::round(y));
#ifndef NDEBUG
if (
p.first < 0 &&
p.first >=
m.size1())
if (
p.second < 0 &&
p.second >=
m.size2())
#endif
}
const auto &intensity =
iI[
i];
}
}
return 1. /
m(idx.first, idx.second);
}
BoundaryOp(boost::shared_ptr<MatrixDouble> flux_ptr,
double &glob_flux)
auto t_flux = getFTensor1FromMat<3>(*
fluxPtr);
for (auto gg = 0; gg != nb_gauss_pts; ++gg) {
++t_flux;
++t_w;
};
}
private:
boost::shared_ptr<MatrixDouble>
fluxPtr;
};
char images_array_files[255] = "out_arrays.txt";
images_array_files, 255, PETSC_NULL);
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;
auto it = data.begin();
if (it == data.end()) {
}
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Number of images " << *it;
it++;
for (; it != data.end();) {
<<
"Read data set " <<
rZ.back() <<
" size " <<
r <<
" by " <<
c;
for (
auto rit =
m.begin1(); rit !=
m.end1(); ++rit) {
for (auto cit = rit.begin(); cit != rit.end(); ++cit) {
*cit = *(it++);
}
}
}
};
structure_data(read_images());
int window_shift = 0;
PETSC_NULL);
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Wave number " <<
k;
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Window shift " << window_shift;
"Expected even number of images");
if (!d1_sg_data)
"Wrong Savitzky Golay order");
if (!d1_sg_data_window)
"Wrong Savitzky Golay window");
}
auto get_bounding_box = [&]() {
MOAB_THROW(moab.get_entities_by_type(0, MBVERTEX, verts));
ParallelComm *pcomm =
CHKERR pcomm->filter_pstatus(verts, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &verts_part);
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);
}
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]};
};
}
1);
int base_order = 1;
PETSC_NULL);
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Base order " << base_order;
}
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
pipeline_mng->setBoundaryRhsIntegrationRule(rule_vol);
auto flux_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
pipeline_mng->getOpBoundaryRhsPipeline().push_back(
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;
}
pipeline_mng->getDomainRhsFE().reset();
pipeline_mng->getBoundaryRhsFE().reset();
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(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
pipeline_mng->getOpDomainLhsPipeline().push_back(
auto unity = []() { return 1; };
pipeline_mng->getOpDomainLhsPipeline().push_back(
new OpHdivU(
"S",
"PHI", unity,
true));
pipeline_mng->getOpDomainRhsPipeline().push_back(
}
CHKERR KSPSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(iD, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(iD, INSERT_VALUES, SCATTER_FORWARD);
double i_lambda_flux;
<< "iD flux " << std::scientific << i_lambda_flux;
}
auto post_proc_fe = boost::make_shared<PostProcEle>(
mField);
auto jac_ptr = boost::make_shared<MatrixDouble>();
post_proc_fe->getOpPtrVector().push_back(
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(
post_proc_fe->getOpPtrVector().push_back(
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{}
)
);
CHKERR post_proc_fe->writeFile(
"out_" + boost::lexical_cast<std::string>(
i) +
".h5m");
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
}
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#define CATCH_ERRORS
Catch errors.
#define MOAB_THROW(err)
Check error code of MoAB function and throw MoFEM exception.
@ L2
field with C-1 continuity
@ HCURL
field with continuous tangents
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
auto createDMVector(DM dm)
Get smart vector from DM.
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.
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
FTensor::Index< 'i', SPACE_DIM > i
const double c
speed of light (cm/ns)
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'k', 3 > k
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMixDivTimesScalar< 2 > OpHdivU
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
const int * d1_normalisation[]
constexpr std::array< int, 3 > d1_savitzky_golay_w3_p2
constexpr std::array< int, 10 > d1_normalisation_p2
constexpr std::array< int, 9 > d1_savitzky_golay_w9_p4
constexpr std::array< int, 7 > d1_savitzky_golay_w7_p4
constexpr std::array< int, 5 > d1_savitzky_golay_w5_p2
constexpr std::array< int, 10 > d1_normalisation_p4
static int window_savitzky_golay
static int order_savitzky_golay
const int ** d1_savitzky_golay[]
constexpr std::array< int, 7 > d1_savitzky_golay_w7_p2
constexpr std::array< int, 5 > d1_savitzky_golay_w5_p4
const int * d1_savitzky_golay_p2[]
const int * d1_savitzky_golay_p4[]
constexpr std::array< int, 9 > d1_savitzky_golay_w9_p2
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
OpBaseImpl< PETSC, EdgeEleOp > OpBase
FTensor::Index< 'i', 3 > i
FTensor::Index< 'm', 3 > m
boost::shared_ptr< MatrixDouble > fluxPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
BoundaryOp(boost::shared_ptr< MatrixDouble > flux_ptr, double &glob_flux)
MoFEMErrorCode readMesh()
[Run problem]
static std::array< double, LAST_BB > aveMaxMin
static std::vector< MatrixInt > iI
static double rhsSource(const double x, const double y, const double)
MoFEMErrorCode solveSystem()
[Solve]
static std::pair< int, int > getCoordsInImage(double x, double y)
MoFEMErrorCode calculateFlux(double &calc_flux)
[Set up problem]
MoFEMErrorCode runProblem()
[Run problem]
MoFEMErrorCode assembleSystemFlux()
MoFEM::Interface & mField
MoFEMErrorCode assembleSystemIntensity()
[Calculate flux on boundary]
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
static double lhsFlux(const double x, const double y, const double)
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode outputResults()
[Solve]
static int savitzkyGolayNormalisation
static std::vector< double > rZ
static const int * savitzkyGolayWeights
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
default operator for EDGE element
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.
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()
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 3, SPACE_DIM > OpHdivHdiv
[Linear elastic problem]