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;
};
};
}
moab::Interface::UNION);
double l2;
for (auto e : edges) {
double l2e = 0;
for (
int j = 0;
j != 3; ++
j) {
l2e += pow(coords(0,
j) - coords(1,
j), 2);
}
l2 = std::max(l2, l2e);
}
}
}
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 dofs_ptr = problem_ptr->getNumeredRowDofsPtr();
const auto uid = DofEntity::getUniqueIdCalculate(
0, FieldEntity::getLocalUniqueIdCalculate(bit_number,
v));
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",
auto solver = pipeline_mng->createKSP();
CHKERR KSPSetFromOptions(solver);
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[]) {
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
}
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);
}
}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#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.
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.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
const double v
phase velocity of light in medium (cm/ns)
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
FormsIntegrators< DomainEleOp >::Assembly< A >::BiLinearForm< I >::OpGradGrad< 1, FIELD_DIM, SPACE_DIM > OpGradGrad
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
boost::shared_ptr< MatrixDouble > uGradPtr
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
OpRhs(boost::shared_ptr< MatrixDouble > u_grad_ptr)
MoFEMErrorCode boundaryCondition()
[Set up problem]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEMErrorCode readMesh()
[Run problem]
MoFEMErrorCode setIntegrationRules()
[Set up problem]
MoFEMErrorCode checkResults()
[Postprocess results]
MoFEMErrorCode solveSystem()
[Solve]
MoFEMErrorCode createCommonData()
[Set up problem]
MoFEMErrorCode runProblem()
[Run problem]
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
[Run problem]
MoFEMErrorCode outputResults()
[Solve]
virtual FieldBitNumber get_field_bit_number(const std::string name) const =0
get field bit number
virtual moab::Interface & get_moab()=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.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Problem manager is used to build and partition problems.
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.