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>(
152 std::array<double, 3> error = {0, 0,
155 for (
int gg = 0; gg != nb_integration_pts; ++gg) {
157 const double alpha = t_w * volume;
161 double diff = t_val -
fun(t_coords(0), t_coords(1), t_coords(2));
162 error[0] += alpha * pow(diff, 2);
163 auto t_diff_grad =
diff_fun(t_coords(0), t_coords(1), t_coords(2));
164 t_diff_grad(
i) -= t_grad_val(
i);
166 error[1] += alpha * t_diff_grad(
i) *
170 MOFEM_LOG(
"SELF", Sev::noisy) <<
"t_hessian_val " << t_hessian_push;
173 if (std::abs(t_hessian_val(0, 1) - t_hessian_val(1, 0)) >
174 std::numeric_limits<float>::epsilon()) {
175 MOFEM_LOG(
"SELF", Sev::error) <<
"t_hessian_val " << t_hessian_val;
177 "Hessian should be symmetric");
180 auto t_diff_hessian =
diff2_fun(t_coords(0), t_coords(1), t_coords(2));
181 t_diff_hessian(
i,
j) -= t_hessian_val(
i,
j);
182 error[2] = t_diff_hessian(
i,
j) * t_diff_hessian(
i,
j);
194 std::array<int, 3> index = {0, 1, 2};
232 constexpr int order = 4;
244 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
265 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
268 CHKERR KSPSetFromOptions(solver);
276 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
277 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
290 auto rule = [](
int,
int,
int p) ->
int {
return 2 *
p; };
295 auto common_data_ptr = boost::make_shared<CommonData>();
298 common_data_ptr->approxVals = boost::make_shared<VectorDouble>();
299 common_data_ptr->approxGradVals = boost::make_shared<MatrixDouble>();
300 common_data_ptr->approxHessianVals = boost::make_shared<MatrixDouble>();
303 auto base_mass = boost::make_shared<MatrixDouble>();
304 auto data_l2 = boost::make_shared<EntitiesFieldData>(MBENTITYSET);
305 auto jac_ptr = boost::make_shared<MatrixDouble>();
306 auto det_ptr = boost::make_shared<VectorDouble>();
307 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
308 common_data_ptr->invJacPtr = inv_jac_ptr;
347 common_data_ptr->approxVals));
352 FIELD_NAME, common_data_ptr->approxGradVals));
359 FIELD_NAME, common_data_ptr->approxHessianVals));
368 CHKERR VecZeroEntries(common_data_ptr->L2Vec);
372 CHKERR VecAssemblyBegin(common_data_ptr->L2Vec);
373 CHKERR VecAssemblyEnd(common_data_ptr->L2Vec);
377 constexpr double eps = 1e-8;
380 CHKERR VecGetArrayRead(common_data_ptr->L2Vec, &array);
382 "Error %6.4e Diff Error %6.4e Diff2 Error %6.4e\n",
383 std::sqrt(array[0]), std::sqrt(array[1]), std::sqrt(array[2]));
384 if (std::sqrt(array[0]) >
eps)
386 if (std::sqrt(array[1]) >
eps)
388 "Wrong function first direcative");
389 if (std::sqrt(array[2]) >
eps)
391 "Wrong function second direcative");
393 CHKERR VecRestoreArrayRead(common_data_ptr->L2Vec, &array);
399int main(
int argc,
char *argv[]) {
407 DMType dm_name =
"DMMOFEM";
412 moab::Core mb_instance;
413 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 smartCreateDMVector(DM dm)
Get smart vector from DM.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_vector< 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
auto createSmartVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart 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)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
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.