Calculating geodetic distances using heat method. For details see,
K. Crane, C. Weischedel, M. Wardetzky, Geodesics in heat: A new approach to computing distance based on heat flow, ACM Transactions on Graphics , vol. 32, no. 5, pp. 152:1-152:11, 2013.
Solving two problems in sequence. Show hot to use form integrators and how to implement user data operator.
static char help[] =
"...\n\n";
private:
OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr);
protected:
boost::shared_ptr<MatrixDouble> uGradPtr;
};
};
}
CHKERR mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes);
pinchNodes.insert(nodes[0]);
CHKERR mField.get_moab().get_adjacencies(pinchNodes, 1,
false, edges,
moab::Interface::UNION);
double l2;
for (auto e : edges) {
CHKERR mField.get_moab().get_connectivity(
Range(e, e), nodes,
false);
CHKERR mField.get_moab().get_coords(nodes, &coords(0, 0));
double l2e = 0;
for (
int j = 0;
j != 3; ++
j) {
l2e += pow(coords(0,
j) - coords(1,
j), 2);
}
l2 = std::max(l2, l2e);
}
}
CHKERR simpleInterface->setUp();
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 * p; };
}
}
}
auto grad_u_ptr = boost::make_shared<MatrixDouble>();
auto set_domain_general = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
auto set_domain_lhs_first = [&](auto &pipeline) {
pipeline.push_back(
new OpGradGrad(
"U",
"U", one));
pipeline.push_back(
new OpMass(
"U",
"U",
rho));
};
auto set_domain_rhs_first = [&](auto &pipeline) {};
auto set_domain_lhs_second = [&](auto &pipeline) {
pipeline.push_back(
new OpGradGrad(
"U",
"U", one));
};
auto set_domain_rhs_second = [&](auto &pipeline) {
pipeline.push_back(
new OpRhs(grad_u_ptr));
};
auto solve_first = [&]() {
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
auto dm = simpleInterface->getDM();
auto problem_ptr = mField.get_problem(
simple->getProblemName());
auto bit_number = mField.get_field_bit_number("U");
auto dofs_ptr = problem_ptr->getNumeredRowDofsPtr();
for (
auto v : pinchNodes) {
if (dof != dofs_ptr->end())
VecSetValue(
F, (*dof)->getPetscGlobalDofIdx(), 1, INSERT_VALUES);
}
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
auto solve_second = [&]() {
CHKERR prb_mng->removeDofsOnEntities(
simple->getProblemName(),
"U",
pinchNodes, 0, 1);
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
auto dm = simpleInterface->getDM();
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
};
auto post_proc = [&](const std::string out_name) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto post_proc_fe =
boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(mField);
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
auto u_ptr = boost::make_shared<VectorDouble>();
auto grad_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(),
{{"U", u_ptr}},
{{"GRAD_U", grad_ptr}},
{}, {}));
CHKERR post_proc_fe->writeFile(out_name);
};
CHKERR post_proc(
"out_heat_method_first.h5m");
CHKERR post_proc(
"out_heat_method_second.h5m");
};
}
}
}
int main(
int argc,
char *argv[]) {
const char param_file[] = "param_file.petsc";
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
}
}
if (nb_dofs) {
auto t_grad = getFTensor1FromMat<3>(*
uGradPtr);
auto nb_base_functions = data.
getN().size2();
auto nb_gauss_pts = getGaussPts().size2();
std::array<double, MAX_DOFS_ON_ENTITY> nf;
std::fill(nf.begin(), &nf[nb_dofs], 0);
auto t_w = getFTensor0IntegrationWeight();
for (size_t gg = 0; gg != nb_gauss_pts; ++gg) {
const auto l2 = t_grad(
i) * t_grad(
i);
if (l2 > std::numeric_limits<double>::epsilon())
t_one(
i) = t_grad(
i) / std::sqrt(l2);
else
size_t bb = 0;
for (; bb != nb_dofs; ++bb) {
nf[bb] -= alpha * t_diff_base(
i) * t_one(
i);
++t_diff_base;
}
for (; bb < nb_base_functions; ++bb) {
++t_diff_base;
}
++t_grad;
}
CHKERR VecSetValues<MoFEM::EssentialBcStorage>(getKSPf(), data, &nf[0],
ADD_VALUES);
}
}