24static char help[] =
"...\n\n";
71 OpRhs(boost::shared_ptr<MatrixDouble> u_grad_ptr);
115 moab::Interface::UNION);
117 for (
auto e : edges) {
123 for (
int j = 0;
j != 3; ++
j) {
124 l2e += pow(coords(0,
j) - coords(1,
j), 2);
126 l2 = std::max(l2, l2e);
141 constexpr int order = 1;
154 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
187 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
191 auto set_domain_general = [&](
auto &pipeline) {
192 auto det_ptr = boost::make_shared<VectorDouble>();
193 auto jac_ptr = boost::make_shared<MatrixDouble>();
194 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
201 auto set_domain_lhs_first = [&](
auto &pipeline) {
203 pipeline.push_back(
new OpGradGrad(
"U",
"U", one));
205 pipeline.push_back(
new OpMass(
"U",
"U",
rho));
211 auto set_domain_rhs_first = [&](
auto &pipeline) {};
215 auto set_domain_lhs_second = [&](
auto &pipeline) {
217 pipeline.push_back(
new OpGradGrad(
"U",
"U", one));
221 auto set_domain_rhs_second = [&](
auto &pipeline) {
223 pipeline.push_back(
new OpRhs(grad_u_ptr));
227 auto solve_first = [&]() {
232 auto solver = pipeline_mng->createKSP();
233 CHKERR KSPSetFromOptions(solver);
244 auto dofs_ptr = problem_ptr->getNumeredRowDofsPtr();
249 if (dof != dofs_ptr->end())
250 VecSetValue(
F, (*dof)->getPetscGlobalDofIdx(), 1, INSERT_VALUES);
257 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
258 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
264 auto solve_second = [&]() {
270 CHKERR prb_mng->removeDofsOnEntities(
simple->getProblemName(),
"U",
275 auto solver = pipeline_mng->createKSP();
276 CHKERR KSPSetFromOptions(solver);
284 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
285 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
291 auto post_proc = [&](
const std::string out_name) {
296 auto det_ptr = boost::make_shared<VectorDouble>();
297 auto jac_ptr = boost::make_shared<MatrixDouble>();
298 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
302 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
304 post_proc_fe->getOpPtrVector().push_back(
306 post_proc_fe->getOpPtrVector().push_back(
308 post_proc_fe->getOpPtrVector().push_back(
311 auto u_ptr = boost::make_shared<VectorDouble>();
312 auto grad_ptr = boost::make_shared<MatrixDouble>();
314 post_proc_fe->getOpPtrVector().push_back(
316 post_proc_fe->getOpPtrVector().push_back(
321 post_proc_fe->getOpPtrVector().push_back(
new OpPPMap(
322 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
326 {{
"GRAD_U", grad_ptr}},
332 CHKERR post_proc_fe->writeFile(out_name);
344 CHKERR post_proc(
"out_heat_method_first.h5m");
355 CHKERR post_proc(
"out_heat_method_second.h5m");
382int main(
int argc,
char *argv[]) {
389 auto core_log = logging::core::get();
398 DMType dm_name =
"DMMOFEM";
403 moab::Core mb_instance;
404 moab::Interface &moab = mb_instance;
434 auto t_grad = getFTensor1FromMat<3>(*uGradPtr);
436 auto nb_base_functions = data.
getN().size2();
437 auto nb_gauss_pts = getGaussPts().size2();
438 std::array<double, MAX_DOFS_ON_ENTITY> nf;
439 std::fill(nf.begin(), &nf[nb_dofs], 0);
442 auto t_w = getFTensor0IntegrationWeight();
443 auto a = getMeasure();
445 for (
size_t gg = 0; gg != nb_gauss_pts; ++gg) {
446 double alpha = t_w *
a;
448 const auto l2 = t_grad(
i) * t_grad(
i);
450 if (l2 > std::numeric_limits<double>::epsilon())
451 t_one(
i) = t_grad(
i) / std::sqrt(l2);
456 for (; bb != nb_dofs; ++bb) {
457 nf[bb] -= alpha * t_diff_base(
i) * t_one(
i);
461 for (; bb < nb_base_functions; ++bb) {
468 CHKERR VecSetValues<MoFEM::EssentialBcStorage>(getKSPf(), data, &nf[0],
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.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#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.
implementation of Data Operators for Forces and Sources
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpCalculateHOJacForFaceImpl< 3 > OpCalculateHOJacForFaceEmbeddedIn3DSpace
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()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
MoFEMErrorCode createCommonData()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
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.
static UId getUniqueIdCalculate(const DofIdx dof, UId ent_uid)
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.
UId getLocalUniqueIdCalculate()
Get the Local Unique Id Calculate object.
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.
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()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
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 addDomainField(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
MoFEMErrorCode getOptions()
get options
MoFEMErrorCode getDM(DM *dm)
Get DM.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.