static char help[] =
"...\n\n";
};
constexpr double a2 =
a *
a;
constexpr double a4 =
a2 *
a2;
constexpr double A = 6371220;
auto res_J = [](
const double x,
const double y,
const double z) {
const double res = (x * x + y * y + z * z -
a2);
return res;
};
auto res_J_dx = [](
const double x,
const double y,
const double z) {
const double res =
res_J(x, y, z);
res * (2 * z)};
};
auto lhs_J_dx2 = [](
const double x,
const double y,
const double z) {
const double res =
res_J(x, y, z);
(res * 2 + (4 * x * x)),
(4 * y * x),
(4 * z * x),
(4 * x * y),
(2 * res + (4 * y * y)),
(4 * z * y),
(4 * x * z),
(4 * y * z),
(2 * res + (4 * z * z))};
};
boost::shared_ptr<MatrixDouble> dot_x_ptr)
auto t_w = getFTensor0IntegrationWeight();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*
xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*
xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
t_P(
i,
j) = t_n(
i) * t_n(
j);
auto t_J_res =
res_J_dx(t_x(0), t_x(1), t_x(2));
const double alpha = t_w;
auto t_nf = getFTensor1FromArray<3, 3>(
locF);
double l = std::sqrt(t_normal(
i) * t_normal(
i));
alpha *
l * ((t_P(
i,
k) * t_J_res(
k) + t_Q(
i,
k) * t_dot_x(
k)));
int rr = 0;
for (; rr !=
nbRows / 3; ++rr) {
t_nf(
j) += t_row_base * t_res(
j);
++t_row_base;
++t_nf;
}
++t_row_base;
}
++t_w;
++t_x;
++t_dot_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble>
xPtr;
boost::shared_ptr<MatrixDouble>
xDotPtr;
};
boost::shared_ptr<MatrixDouble> dot_x_ptr)
this->sYmm = false;
}
auto t_w = getFTensor0IntegrationWeight();
auto t_x0 = getFTensor1CoordsAtGaussPts();
auto t_x = getFTensor1FromMat<3>(*
xPtr);
auto t_dot_x = getFTensor1FromMat<3>(*
xDotPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
auto get_t_mat = [&](const int rr) {
};
const double ts_a = getTSa();
t_P(
i,
j) = t_n(
i) * t_n(
j);
auto t_J_lhs =
lhs_J_dx2(t_x(0), t_x(1), t_x(2));
double l = std::sqrt(t_normal(
i) * t_normal(
i));
const double alpha = t_w;
(
alpha *
l) * (t_P(
i,
k) * t_J_lhs(
k,
j) + t_Q(
i,
j) * ts_a);
int rr = 0;
for (; rr !=
nbRows / 3; rr++) {
auto t_mat = get_t_mat(3 * rr);
for (
int cc = 0; cc !=
nbCols / 3; cc++) {
const double rc = t_row_base * t_col_base;
t_mat(
i,
j) += rc * t_lhs(
i,
j);
++t_col_base;
++t_mat;
}
++t_row_base;
}
++t_row_base;
++t_w;
++t_x;
++t_x0;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble>
xPtr;
boost::shared_ptr<MatrixDouble>
xDotPtr;
};
}
auto t_x = getFTensor1FromMat<3>(*
xPtr);
auto t_normal = getFTensor1NormalsAtGaussPts();
double error = 0;
for (int gg = 0; gg != nb_integration_pts; gg++) {
double l = std::sqrt(t_normal(
i) * t_normal(
i));
error += t_w *
l * std::abs((t_x(
i) * t_x(
i) -
A *
A));
++t_w;
++t_x;
++t_normal;
}
}
private:
boost::shared_ptr<MatrixDouble>
xPtr;
};
private:
};
}
}
}
}
};
auto x_ptr = boost::make_shared<MatrixDouble>();
auto dot_x_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
auto def_ops = [&](auto &pipeline) {
pipeline.push_back(
pipeline.push_back(
};
new OpRhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
new OpLhs(
"HO_POSITIONS", x_ptr, dot_x_ptr));
}
ts = pipeline_mng->createTSIM();
double ftime = 1;
CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
SCATTER_FORWARD);
}
auto x_ptr = boost::make_shared<MatrixDouble>();
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(
new OpPPMap(post_proc_fe->getPostProcMesh(),
post_proc_fe->getMapGaussPts(),
{},
{{"HO_POSITIONS", x_ptr}},
{}, {}
)
);
CHKERR post_proc_fe->writeFile(
"out_approx.h5m");
auto error_fe = boost::make_shared<DomainEle>(
mField);
error_fe->getOpPtrVector().push_back(
error_fe->getOpPtrVector().push_back(
error_fe->getOpPtrVector().push_back(
new OpError(
"HO_POSITIONS", x_ptr));
error_fe->preProcessHook = [&]() {
MOFEM_LOG(
"EXAMPLE", Sev::inform) <<
"Create vec ";
};
error_fe->postProcessHook = [&]() {
double error2;
<<
"Error " << std::sqrt(error2 / (4 * M_PI *
A *
A));
};
"PARALLEL=WRITE_PART");
else
}
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;
}
}
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.
#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< 'i', SPACE_DIM > i
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
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]
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
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.
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
boost::shared_ptr< MatrixDouble > xPtr
boost::shared_ptr< MatrixDouble > xDotPtr
MoFEMErrorCode iNtegrate(EntitiesFieldData::EntData &row_data)
Class dedicated to integrate operator.