23 using namespace MoFEM;
25 static char help[] =
"...\n\n";
27 #include <BasicFiniteElements.hpp>
46 constexpr
double a = 1;
47 constexpr
double a2 =
a *
a;
50 constexpr
double A = 6371220;
56 auto res_J = [](
const double x,
const double y,
const double z) {
57 const double res = (x * x + y * y + z * z -
a2);
61 auto res_J_dx = [](
const double x,
const double y,
const double z) {
62 const double res =
res_J(x, y, z);
67 auto lhs_J_dx2 = [](
const double x,
const double y,
const double z) {
68 const double res =
res_J(x, y, z);
71 (res * 2 + (4 * x * x)),
76 (2 * res + (4 * y * y)),
81 (2 * res + (4 * z * z))};
86 OpRhs(
const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
87 boost::shared_ptr<MatrixDouble> dot_x_ptr)
89 xPtr(x_ptr), xDotPtr(dot_x_ptr) {}
94 auto t_w = getFTensor0IntegrationWeight();
97 auto t_x0 = getFTensor1CoordsAtGaussPts();
98 auto t_x = getFTensor1FromMat<3>(*xPtr);
99 auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
100 auto t_normal = getFTensor1NormalsAtGaussPts();
104 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
110 t_P(
i,
j) = t_n(
i) * t_n(
j);
113 auto t_J_res =
res_J_dx(t_x(0), t_x(1), t_x(2));
115 const double alpha = t_w;
116 auto t_nf = getFTensor1FromArray<3, 3>(locF);
117 double l = std::sqrt(t_normal(
i) * t_normal(
i));
121 alpha *
l * ((t_P(
i,
k) * t_J_res(
k) + t_Q(
i,
k) * t_dot_x(
k)));
124 for (; rr != nbRows / 3; ++rr) {
126 t_nf(
j) += t_row_base * t_res(
j);
131 for (; rr < nbRowBaseFunctions; ++rr) {
146 boost::shared_ptr<MatrixDouble>
xPtr;
152 OpLhs(
const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr,
153 boost::shared_ptr<MatrixDouble> dot_x_ptr)
156 xPtr(x_ptr), xDotPtr(dot_x_ptr) {
164 auto t_w = getFTensor0IntegrationWeight();
167 auto t_x0 = getFTensor1CoordsAtGaussPts();
168 auto t_x = getFTensor1FromMat<3>(*xPtr);
169 auto t_dot_x = getFTensor1FromMat<3>(*xDotPtr);
170 auto t_normal = getFTensor1NormalsAtGaussPts();
172 auto get_t_mat = [&](
const int rr) {
174 &locMat(rr + 0, 0), &locMat(rr + 0, 1), &locMat(rr + 0, 2),
176 &locMat(rr + 1, 0), &locMat(rr + 1, 1), &locMat(rr + 1, 2),
178 &locMat(rr + 2, 0), &locMat(rr + 2, 1), &locMat(rr + 2, 2)};
182 const double ts_a = getTSa();
184 for (
int gg = 0; gg != nbIntegrationPts; gg++) {
190 t_P(
i,
j) = t_n(
i) * t_n(
j);
193 auto t_J_lhs =
lhs_J_dx2(t_x(0), t_x(1), t_x(2));
194 double l = std::sqrt(t_normal(
i) * t_normal(
i));
196 const double alpha = t_w;
199 (
alpha *
l) * (t_P(
i,
k) * t_J_lhs(
k,
j) + t_Q(
i,
j) * ts_a);
202 for (; rr != nbRows / 3; rr++) {
205 auto t_mat = get_t_mat(3 * rr);
207 for (
int cc = 0; cc != nbCols / 3; cc++) {
209 const double rc = t_row_base * t_col_base;
210 t_mat(
i,
j) += rc * t_lhs(
i,
j);
218 for (; rr < nbRowBaseFunctions; ++rr)
231 boost::shared_ptr<MatrixDouble>
xPtr;
237 OpError(
const std::string field_name, boost::shared_ptr<MatrixDouble> x_ptr)
241 std::fill(&doEntities[MBEDGE], &doEntities[MBMAXTYPE],
false);
248 auto t_w = getFTensor0IntegrationWeight();
249 auto t_x = getFTensor1FromMat<3>(*xPtr);
250 auto t_normal = getFTensor1NormalsAtGaussPts();
251 auto nb_integration_pts = getGaussPts().size2();
255 for (
int gg = 0; gg != nb_integration_pts; gg++) {
257 double l = std::sqrt(t_normal(
i) * t_normal(
i));
258 error += t_w *
l * std::abs((t_x(
i) * t_x(
i) -
A *
A));
265 CHKERR VecSetValue(errorVec, 0, error, ADD_VALUES);
273 boost::shared_ptr<MatrixDouble>
xPtr;
351 auto x_ptr = boost::make_shared<MatrixDouble>();
352 auto dot_x_ptr = boost::make_shared<MatrixDouble>();
353 auto det_ptr = boost::make_shared<VectorDouble>();
354 auto jac_ptr = boost::make_shared<MatrixDouble>();
355 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
357 auto def_ops = [&](
auto &pipeline) {
368 new OpRhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
370 new OpLhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
382 CHKERR mField.loop_dofs(
"HO_POSITIONS", ent_method_material);
387 auto dm =
simple->getDM();
389 ts = pipeline_mng->createTSIM();
392 CHKERR TSSetMaxSteps(ts, 1);
393 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
399 CHKERR TSSetFromOptions(ts);
402 CHKERR TSGetTime(ts, &ftime);
414 auto dm =
simple->getDM();
416 auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
417 post_proc_fe->generateReferenceElementMesh();
418 post_proc_fe->addFieldValuesPostProc(
"HO_POSITIONS");
420 CHKERR post_proc_fe->writeFile(
"out_approx.h5m");
422 auto error_fe = boost::make_shared<DomainEle>(mField);
424 auto x_ptr = boost::make_shared<MatrixDouble>();
425 auto det_ptr = boost::make_shared<VectorDouble>();
426 auto jac_ptr = boost::make_shared<MatrixDouble>();
427 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
429 error_fe->getOpPtrVector().push_back(
431 error_fe->getOpPtrVector().push_back(
433 error_fe->getOpPtrVector().push_back(
new OpError(
"HO_POSITIONS", x_ptr));
435 error_fe->preProcessHook = [&]() {
437 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Create vec ";
439 mField.get_comm(), (!mField.get_comm_rank()) ? 1 : 0, 1);
444 error_fe->postProcessHook = [&]() {
451 <<
"Error " << std::sqrt(error2 / (4 * M_PI *
A *
A));
460 if (mField.get_comm_size() > 1)
461 CHKERR mField.get_moab().write_file(
"out_ho_mesh.h5m",
"MOAB",
462 "PARALLEL=WRITE_PART");
464 CHKERR mField.get_moab().write_file(
"out_ho_mesh.h5m");
469 int main(
int argc,
char *argv[]) {
475 auto core_log = logging::core::get();
484 DMType dm_name =
"DMMOFEM";
void simple(double P1[], double P2[], double P3[], double c[], const int N)
int main(int argc, char *argv[])
[Postprocess results]
EntitiesFieldData::EntData EntData
FTensor::Index< 'j', 3 > j
FTensor::Index< 'i', 3 > i
FTensor::Index< 'k', 3 > k
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
Tensor1< T, Tensor_Dim > normalize()
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#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.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
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.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'l', 3 > l
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
auto createSmartVectorMPI
Create MPI Vector.
DeprecatedCoreInterface Interface
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEM::Interface & mField
MoFEMErrorCode setupProblem()
[Read mesh]
ApproxSphere(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
[Run programme]
MoFEMErrorCode setOPs()
[Set up problem]
MoFEMErrorCode getOptions()
[Run programme]
MoFEMErrorCode outputResults()
[Solve]
PipelineManager::FaceEle DomainEle
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::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
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.
Approximate field valuse for given petsc vector.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Calculate normals at Gauss points of triangle element.
PipelineManager interface.
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
OpError(const std::string field_name, boost::shared_ptr< MatrixDouble > x_ptr)
static SmartPetscObj< Vec > errorVec
boost::shared_ptr< MatrixDouble > xPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data, EntitiesFieldData::EntData &col_data)
Integrate grad-grad operator.
boost::shared_ptr< MatrixDouble > xPtr
boost::shared_ptr< MatrixDouble > xDotPtr
OpLhs(const std::string field_name, boost::shared_ptr< MatrixDouble > x_ptr, boost::shared_ptr< MatrixDouble > dot_x_ptr)
boost::shared_ptr< MatrixDouble > xPtr
boost::shared_ptr< MatrixDouble > xDotPtr
OpRhs(const std::string field_name, boost::shared_ptr< MatrixDouble > x_ptr, boost::shared_ptr< MatrixDouble > dot_x_ptr)
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Class dedicated to integrate operator.