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>;
PETSC>::LinearForm<GAUSS>::OpSource<1, 1>;
private:
boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
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;
};
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())
if (p.second < 0 && p.second >=
m.size2())
#endif
return p;
}
const auto idx = getCoordsInImage(x, y);
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];
}
const auto idx = getCoordsInImage(x, y);
const auto &
m = iI[focalIndex];
return 1. /
m(idx.first, idx.second);
}
BoundaryOp(boost::shared_ptr<MatrixDouble> flux_ptr, double &glob_flux)
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;
};
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()) {
}
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);
<<
"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;
if ((rZ.size() - 1) % 2)
"Expected even number of images");
focalIndex = (rZ.size() - 1) / 2 + window_shift;
MOFEM_LOG(
"WORLD", Sev::inform) <<
"zR for mid plane " << rZ[focalIndex];
if (!d1_sg_data)
"Wrong Savitzky Golay order");
if (!d1_sg_data_window)
"Wrong Savitzky Golay window");
savitzkyGolayNormalisation =
savitzkyGolayWeights = d1_sg_data_window;
CHKERR assembleSystemIntensity();
}
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]};
};
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
aveMaxMin = get_bounding_box();
<< "Centre " << aveMaxMin[CENTER_X] << " " << aveMaxMin[CENTER_Y];
<< "Max " << aveMaxMin[MAX_X] << " " << aveMaxMin[MAX_Y];
<< "Min " << aveMaxMin[MIN_X] << " " << aveMaxMin[MIN_Y];
}
1);
int base_order = 1;
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();
CHKERR simpleInterface->addFieldToEmptyFieldBlocks(
"PHI",
"PHI");
}
auto flux_ptr = boost::make_shared<MatrixDouble>();
new BoundaryOp(flux_ptr, calc_flux));
calc_flux = 0;
double global_flux_assembeld = 0;
MPI_Allreduce(&calc_flux, &global_flux_assembeld, 1, MPI_DOUBLE, MPI_SUM,
calc_flux = global_flux_assembeld;
}
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto unity = []() { return 1; };
new OpHdivU(
"S",
"PHI", unity,
true));
}
auto dm = simpleInterface->getDM();
CHKERR KSPSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(iD, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(iD, INSERT_VALUES, SCATTER_FORWARD);
double i_lambda_flux;
CHKERR calculateFlux(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(),
)
);
CHKERR post_proc_fe->writeFile(
"out_" + boost::lexical_cast<std::string>(
i) +
".h5m");
}
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
try {
DMType dm_name = "DMMOFEM";
}
}