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";
#include <BasicFiniteElements.hpp>
private:
Simple *simpleInterface;
Range pinchNodes;
OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr);
protected:
boost::shared_ptr<MatrixDouble> uGradPtr;
};
};
}
CHKERR mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
Range nodes;
CHKERR mField.get_moab().get_entities_by_type(0, MBVERTEX, nodes);
pinchNodes.insert(nodes[0]);
Range edges;
CHKERR mField.get_moab().get_adjacencies(pinchNodes, 1,
false, edges,
moab::Interface::UNION);
double l2;
for (auto e : edges) {
Range nodes;
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; };
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule);
CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule);
}
}
}
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
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>();
pipeline.push_back(new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
pipeline.push_back(new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
};
auto set_domain_lhs_first = [&](auto &pipeline) {
auto one = [](double, double, double) -> double { return 1.; };
pipeline.push_back(
new OpGradGrad(
"U",
"U", one));
auto rho = [](double, double, double) ->
double {
return 1. /
dt; };
pipeline.push_back(
new OpMass(
"U",
"U",
rho));
};
auto set_domain_rhs_first = [&](auto &pipeline) {};
auto set_domain_lhs_second = [&](auto &pipeline) {
auto one = [](double, double, double) { return 1.; };
pipeline.push_back(
new OpGradGrad(
"U",
"U", one));
};
auto set_domain_rhs_second = [&](auto &pipeline) {
pipeline.push_back(new OpCalculateScalarFieldGradient<3>("U", grad_u_ptr));
pipeline.push_back(
new OpRhs(grad_u_ptr));
};
auto solve_first = [&]() {
auto simple = mField.getInterface<Simple>();
auto pipeline_mng = mField.getInterface<PipelineManager>();
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) {
auto dof = dofs_ptr->get<Unique_mi_tag>().find(uid);
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 = [&]() {
auto simple = mField.getInterface<Simple>();
auto prb_mng = mField.getInterface<ProblemsManager>();
CHKERR prb_mng->removeDofsOnEntities(
simple->getProblemName(),
"U",
pinchNodes, 0, 1);
auto *pipeline_mng = mField.getInterface<PipelineManager>();
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) {
PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
auto tmp_lhs_fe = pipeline_mng->getDomainLhsFE();
auto tmp_rhs_fe = pipeline_mng->getDomainRhsFE();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
pipeline_mng->getDomainLhsFE().reset();
pipeline_mng->getDomainRhsFE().reset();
auto post_proc_fe =
boost::make_shared<PostProcFaceOnRefinedMeshFor2D>(mField);
post_proc_fe->generateReferenceElementMesh();
post_proc_fe->getOpPtrVector().push_back(
post_proc_fe->getOpPtrVector().push_back(
new OpInvertMatrix<3>(jac_ptr, det_ptr, inv_jac_ptr));
post_proc_fe->getOpPtrVector().push_back(
new OpSetInvJacH1ForFaceEmbeddedIn3DSpace(inv_jac_ptr));
post_proc_fe->addFieldValuesPostProc("U");
post_proc_fe->addFieldValuesGradientPostProc("U");
pipeline_mng->getDomainRhsFE() = post_proc_fe;
CHKERR pipeline_mng->loopFiniteElements();
CHKERR post_proc_fe->writeFile(out_name);
tmp_lhs_fe = pipeline_mng->getDomainLhsFE() = tmp_lhs_fe;
tmp_rhs_fe = pipeline_mng->getDomainRhsFE() = tmp_rhs_fe;
};
set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
set_domain_lhs_first(pipeline_mng->getOpDomainLhsPipeline());
set_domain_rhs_first(pipeline_mng->getOpDomainRhsPipeline());
CHKERR post_proc(
"out_heat_method_first.h5m");
pipeline_mng->getOpDomainLhsPipeline().clear();
pipeline_mng->getOpDomainRhsPipeline().clear();
set_domain_general(pipeline_mng->getOpDomainLhsPipeline());
set_domain_general(pipeline_mng->getOpDomainRhsPipeline());
set_domain_lhs_second(pipeline_mng->getOpDomainLhsPipeline());
set_domain_rhs_second(pipeline_mng->getOpDomainRhsPipeline());
CHKERR post_proc(
"out_heat_method_second.h5m");
};
}
}
}
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
try {
DMType dm_name = "DMMOFEM";
}
}
auto nb_dofs = data.getIndices().size();
if (nb_dofs) {
auto t_grad = getFTensor1FromMat<3>(*
uGradPtr);
auto nb_base_functions = data.getN().size2();
std::array<double, MAX_DOFS_ON_ENTITY> nf;
std::fill(nf.begin(), &nf[nb_dofs], 0);
auto t_diff_base = data.getFTensor1DiffN<3>();
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);
}
}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpMass
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto smartCreateDMVector
Get smart vector from DM.
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.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
int main(int argc, char *argv[])
[Check results]
EntitiesFieldData::EntData EntData
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 3 > OpGradGrad
double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'j', 3 > j
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
boost::shared_ptr< MatrixDouble > uGradPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpRhs(boost::shared_ptr< MatrixDouble > u_grad_ptr)
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode setIntegrationRules()
[Set up problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Create common data]
MoFEMErrorCode runProblem()
[Create common data]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode outputResults()
[Solve]
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.
static UId getUniqueIdCalculate(const DofIdx dof, UId ent_uid)
Data on single entity (This is passed as argument to DataOperator::doWork)
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator