static char help[] =
"...\n\n";
};
auto fun = [](
const double x,
const double y,
const double z) {
return x * x + y * y + x * y + pow(x, 3) + pow(y, 3) + pow(x, 4) + pow(y, 4);
};
auto diff_fun = [](
const double x,
const double y,
const double z) {
2 * x + y + 3 * pow(x, 2) + 4 * pow(x, 3),
2 * y + x + 3 * pow(y, 2) + 4 * pow(y, 3)};
};
auto diff2_fun = [](
const double x,
const double y,
const double z) {
2 + 6 * x + 12 * pow(x, 2), 1.,
1., 2 + 6 * y + 12 * pow(y, 2)};
};
private:
};
};
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
}
auto t_val = getFTensor0FromVec(*(
auto t_grad_val = getFTensor1FromMat<SPACE_DIM>(
->approxGradVals));
auto t_hessian_val = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
auto t_inv_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
std::array<double, 3> error = {0, 0,
0};
for (int gg = 0; gg != nb_integration_pts; ++gg) {
const double alpha = t_w * volume;
double diff = t_val -
fun(t_coords(0), t_coords(1), t_coords(2));
error[0] +=
alpha * pow(diff, 2);
auto t_diff_grad =
diff_fun(t_coords(0), t_coords(1), t_coords(2));
t_diff_grad(
i) -= t_grad_val(
i);
error[1] +=
alpha * t_diff_grad(
i) *
MOFEM_LOG(
"SELF", Sev::noisy) <<
"t_hessian_val " << t_hessian_push;
if (std::abs(t_hessian_val(0, 1) - t_hessian_val(1, 0)) >
std::numeric_limits<float>::epsilon()) {
MOFEM_LOG(
"SELF", Sev::error) <<
"t_hessian_val " << t_hessian_val;
"Hessian should be symmetric");
}
auto t_diff_hessian =
diff2_fun(t_coords(0), t_coords(1), t_coords(2));
t_diff_hessian(
i,
j) -= t_hessian_val(
i,
j);
error[2] = t_diff_hessian(
i,
j) * t_diff_hessian(
i,
j);
++t_w;
++t_val;
++t_grad_val;
++t_hessian_val;
++t_inv_jac;
++t_coords;
}
std::array<int, 3>
index = {0, 1, 2};
ADD_VALUES);
}
};
}
}
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
}
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
CHKERR KSPSetFromOptions(solver);
CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
rule);
auto common_data_ptr = boost::make_shared<CommonData>();
common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
common_data_ptr->approxGradVals = boost::make_shared<MatrixDouble>();
common_data_ptr->approxHessianVals = boost::make_shared<MatrixDouble>();
auto base_mass = boost::make_shared<MatrixDouble>();
auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto det_ptr = boost::make_shared<VectorDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
common_data_ptr->invJacPtr = inv_jac_ptr;
base_mass, data_l2,
base_mass, data_l2,
base_mass, data_l2,
common_data_ptr->approxVals));
CHKERR VecZeroEntries(common_data_ptr->L2Vec);
CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
constexpr double eps = 1e-8;
const double *array;
CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
"Error %6.4e Diff Error %6.4e Diff2 Error %6.4e\n",
std::sqrt(array[0]), std::sqrt(array[1]), std::sqrt(array[2]));
if (std::sqrt(array[0]) >
eps)
if (std::sqrt(array[1]) >
eps)
"Wrong function first direcative");
if (std::sqrt(array[2]) >
eps)
"Wrong function second direcative");
CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
}
int main(
int argc,
char *argv[]) {
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
}
}
#define MOFEM_LOG_C(channel, severity, format,...)
constexpr char FIELD_NAME[]
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
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
MoFEMErrorCode VecSetValues(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Assemble PETSc vector.
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.
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()
[Check results]
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode readMesh()
[Run programme]
MoFEMErrorCode assembleSystem()
[Push operators to pipeline]
MoFEM::Interface & mField
MoFEMErrorCode runProblem()
[Run programme]
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)
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()
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
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.