14static char help[] =
"...\n\n";
37auto fun = [](
const double x,
const double y,
const double z) {
38 return x * x + y * y + x * y + pow(x, 3) + pow(y, 3) + pow(x, 4) + pow(y, 4);
45auto diff_fun = [](
const double x,
const double y,
const double z) {
47 2 * x + y + 3 * pow(x, 2) + 4 * pow(x, 3),
48 2 * y + x + 3 * pow(y, 2) + 4 * pow(y, 3)};
55auto diff2_fun = [](
const double x,
const double y,
const double z) {
57 2 + 6 * x + 12 * pow(x, 2), 1.,
59 1., 2 + 6 * y + 12 * pow(y, 2)};
117 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
125 const int nb_integration_pts =
getGaussPts().size2();
129 auto t_grad_val = getFTensor1FromMat<SPACE_DIM>(
132 auto t_hessian_val = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
136 auto t_inv_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
150 std::array<double, 3> error = {0, 0,
153 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
155 const double alpha = t_w * volume;
159 double diff = t_val -
fun(t_coords(0), t_coords(1), t_coords(2));
160 error[0] += alpha * pow(diff, 2);
161 auto t_diff_grad =
diff_fun(t_coords(0), t_coords(1), t_coords(2));
162 t_diff_grad(
i) -= t_grad_val(
i);
164 error[1] += alpha * t_diff_grad(
i) *
168 MOFEM_LOG(
"SELF", Sev::noisy) <<
"t_hessian_val " << t_hessian_push;
171 if (std::abs(t_hessian_val(0, 1) - t_hessian_val(1, 0)) >
172 std::numeric_limits<float>::epsilon()) {
173 MOFEM_LOG(
"SELF", Sev::error) <<
"t_hessian_val " << t_hessian_val;
175 "Hessian should be symmetric");
178 auto t_diff_hessian =
diff2_fun(t_coords(0), t_coords(1), t_coords(2));
179 t_diff_hessian(
i,
j) -= t_hessian_val(
i,
j);
180 error[2] = t_diff_hessian(
i,
j) * t_diff_hessian(
i,
j);
192 std::array<int, 3> index = {0, 1, 2};
230 constexpr int order = 4;
242 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
263 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
266 CHKERR KSPSetFromOptions(solver);
274 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
275 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
288 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
293 auto common_data_ptr = boost::make_shared<CommonData>();
296 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
297 common_data_ptr->approxGradVals = boost::make_shared<MatrixDouble>();
298 common_data_ptr->approxHessianVals = boost::make_shared<MatrixDouble>();
301 auto base_mass = boost::make_shared<MatrixDouble>();
302 auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
303 auto jac_ptr = boost::make_shared<MatrixDouble>();
304 auto det_ptr = boost::make_shared<VectorDouble>();
305 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
306 common_data_ptr->invJacPtr = inv_jac_ptr;
345 common_data_ptr->approxVals));
350 FIELD_NAME, common_data_ptr->approxGradVals));
357 FIELD_NAME, common_data_ptr->approxHessianVals));
366 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
370 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
371 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
375 constexpr double eps = 1e-8;
378 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
380 "Error %6.4e Diff Error %6.4e Diff2 Error %6.4e\n",
381 std::sqrt(array[0]), std::sqrt(array[1]), std::sqrt(array[2]));
382 if (std::sqrt(array[0]) >
eps)
384 if (std::sqrt(array[1]) >
eps)
386 "Wrong function first direcative");
387 if (std::sqrt(array[2]) >
eps)
389 "Wrong function second direcative");
391 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
397int main(
int argc,
char *argv[]) {
405 DMType dm_name =
"DMMOFEM";
410 moab::Core mb_instance;
411 moab::Interface &moab = mb_instance;
#define MOFEM_LOG_C(channel, severity, format,...)
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, FIELD_DIM > OpDomainSource
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
@ L2
field with C-1 continuity
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
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.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_deque< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_deque< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
FTensor::Index< 'i', SPACE_DIM > i
auto diff2_fun
Function second derivative.
constexpr char FIELD_NAME[]
auto fun
Function to approximate.
auto diff_fun
Function derivative.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< BASE_DIM, FIELD_DIM > OpDomainMass
OPerator to integrate mass matrix for least square approximation.
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpDomainSource
Operator to integrate the right hand side matrix for the problem.
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
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
Collected data use d by operator to evaluate errors for the test.
boost::shared_ptr< VectorDouble > approxVals
boost::shared_ptr< MatrixDouble > approxGradVals
SmartPetscObj< Vec > L2Vec
boost::shared_ptr< MatrixDouble > invJacPtr
boost::shared_ptr< MatrixDouble > approxHessianVals
Operator to evaluate errors.
OpError(boost::shared_ptr< CommonData > &common_data_ptr)
boost::shared_ptr< CommonData > commonDataPtr
MoFEMErrorCode doWork(int side, EntityType type, EntData &data)
Operator for linear form, usually to calculate values on right hand side.
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
AtomTest(MoFEM::Interface &m_field)
MoFEMErrorCode setupProblem()
MoFEMErrorCode readMesh()
MoFEMErrorCode assembleSystem()
MoFEM::Interface & mField
MoFEMErrorCode runProblem()
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)
default operator for TRI element
auto getFTensor1CoordsAtGaussPts()
Get coordinates at integration points assuming linear geometry.
auto getFTensor0IntegrationWeight()
Get integration weights.
double getMeasure() const
get measure of element
@ OPROW
operator doWork function is executed on FE rows
MatrixDouble & getGaussPts()
matrix of integration (Gauss) points for Volume Element
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Evaluate field gradient values for scalar field, i.e. gradient is tensor rank 1 (vector),...
Get value at integration points for scalar field.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
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.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.