11static char help[] =
"...\n\n";
30constexpr double a = 1;
31constexpr double a2 =
a *
a;
34constexpr double A = 6371220;
40auto res_J = [](
const double x,
const double y,
const double z) {
41 const double res = (x * x + y * y + z * z -
a2);
45auto res_J_dx = [](
const double x,
const double y,
const double z) {
46 const double res =
res_J(x, y, z);
51auto lhs_J_dx2 = [](
const double x,
const double y,
const double z) {
52 const double res =
res_J(x, y, z);
55 (res * 2 + (4 * x * x)),
60 (2 * res + (4 * y * y)),
65 (2 * res + (4 * z * z))};
71 boost::shared_ptr<MatrixDouble> dot_x_ptr)
78 auto t_w = getFTensor0IntegrationWeight();
81 auto t_x0 = getFTensor1CoordsAtGaussPts();
82 auto t_x = getFTensor1FromMat<3>(*
xPtr);
83 auto t_dot_x = getFTensor1FromMat<3>(*
xDotPtr);
84 auto t_normal = getFTensor1NormalsAtGaussPts();
94 t_P(
i,
j) = t_n(
i) * t_n(
j);
97 auto t_J_res =
res_J_dx(t_x(0), t_x(1), t_x(2));
99 const double alpha = t_w;
100 auto t_nf = getFTensor1FromArray<3, 3>(
locF);
101 double l = std::sqrt(t_normal(
i) * t_normal(
i));
105 alpha *
l * ((t_P(
i,
k) * t_J_res(
k) + t_Q(
i,
k) * t_dot_x(
k)));
108 for (; rr !=
nbRows / 3; ++rr) {
110 t_nf(
j) += t_row_base * t_res(
j);
130 boost::shared_ptr<MatrixDouble>
xPtr;
137 boost::shared_ptr<MatrixDouble> dot_x_ptr)
148 auto t_w = getFTensor0IntegrationWeight();
151 auto t_x0 = getFTensor1CoordsAtGaussPts();
152 auto t_x = getFTensor1FromMat<3>(*
xPtr);
153 auto t_dot_x = getFTensor1FromMat<3>(*
xDotPtr);
154 auto t_normal = getFTensor1NormalsAtGaussPts();
156 auto get_t_mat = [&](
const int rr) {
166 const double ts_a = getTSa();
174 t_P(
i,
j) = t_n(
i) * t_n(
j);
177 auto t_J_lhs =
lhs_J_dx2(t_x(0), t_x(1), t_x(2));
178 double l = std::sqrt(t_normal(
i) * t_normal(
i));
180 const double alpha = t_w;
183 (alpha *
l) * (t_P(
i,
k) * t_J_lhs(
k,
j) + t_Q(
i,
j) * ts_a);
186 for (; rr !=
nbRows / 3; rr++) {
189 auto t_mat = get_t_mat(3 * rr);
191 for (
int cc = 0; cc !=
nbCols / 3; cc++) {
193 const double rc = t_row_base * t_col_base;
194 t_mat(
i,
j) += rc * t_lhs(
i,
j);
215 boost::shared_ptr<MatrixDouble>
xPtr;
233 auto t_x = getFTensor1FromMat<3>(*
xPtr);
234 auto t_normal = getFTensor1NormalsAtGaussPts();
239 for (
int gg = 0; gg != nb_integration_pts; gg++) {
241 double l = std::sqrt(t_normal(
i) * t_normal(
i));
242 error += t_w *
l * std::abs((t_x(
i) * t_x(
i) -
A *
A));
257 boost::shared_ptr<MatrixDouble>
xPtr;
335 auto x_ptr = boost::make_shared<MatrixDouble>();
336 auto dot_x_ptr = boost::make_shared<MatrixDouble>();
337 auto det_ptr = boost::make_shared<VectorDouble>();
338 auto jac_ptr = boost::make_shared<MatrixDouble>();
339 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
341 auto def_ops = [&](
auto &pipeline) {
352 new OpRhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
354 new OpLhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
371 auto dm =
simple->getDM();
373 ts = pipeline_mng->createTSIM();
376 CHKERR TSSetMaxSteps(ts, 1);
377 CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
383 CHKERR TSSetFromOptions(ts);
386 CHKERR TSGetTime(ts, &ftime);
397 auto x_ptr = boost::make_shared<MatrixDouble>();
398 auto det_ptr = boost::make_shared<VectorDouble>();
399 auto jac_ptr = boost::make_shared<MatrixDouble>();
400 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
403 auto dm =
simple->getDM();
406 boost::make_shared<PostProcBrokenMeshInMoab<DomainEle>>(
mField);
408 post_proc_fe->getOpPtrVector().push_back(
413 post_proc_fe->getOpPtrVector().push_back(
415 new OpPPMap(post_proc_fe->getPostProcMesh(),
416 post_proc_fe->getMapGaussPts(),
420 {{
"HO_POSITIONS", x_ptr}},
429 CHKERR post_proc_fe->writeFile(
"out_approx.h5m");
431 auto error_fe = boost::make_shared<DomainEle>(
mField);
433 error_fe->getOpPtrVector().push_back(
435 error_fe->getOpPtrVector().push_back(
437 error_fe->getOpPtrVector().push_back(
new OpError(
"HO_POSITIONS", x_ptr));
439 error_fe->preProcessHook = [&]() {
441 MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Create vec ";
448 error_fe->postProcessHook = [&]() {
455 <<
"Error " << std::sqrt(error2 / (4 * M_PI *
A *
A));
466 "PARALLEL=WRITE_PART");
473int main(
int argc,
char *argv[]) {
479 auto core_log = logging::core::get();
488 DMType dm_name =
"DMMOFEM";
493 moab::Core mb_instance;
494 moab::Interface &moab = mb_instance;
void simple(double P1[], double P2[], double P3[], double c[], const int N)
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 ...
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.
auto createDMVector(DM dm)
Get smart vector from DM.
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(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
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 createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
constexpr auto field_name
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]
virtual int get_comm_size() const =0
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =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.
std::array< bool, MBMAXTYPE > doEntities
If true operator is executed for entity.
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.
auto getFTensor0IntegrationWeight()
Get integration weights.
@ OPROW
operator doWork function is executed on FE rows
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.
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
int nbIntegrationPts
number of integration points
MatrixDouble locMat
local entity block matrix
int nbCols
number if dof on column
int nbRowBaseFunctions
number or row base functions
Approximate field values 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.
Post post-proc data at points from hash maps.
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.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
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.