12 using namespace MoFEM;
14 static char help[] =
"...\n\n";
37 auto 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);
45 auto 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)};
55 auto 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)};
97 boost::shared_ptr<MatrixDouble> invJacPtr;
98 boost::shared_ptr<VectorDouble> approxVals;
99 boost::shared_ptr<MatrixDouble> approxGradVals;
100 boost::shared_ptr<MatrixDouble> approxHessianVals;
116 boost::shared_ptr<CommonData> commonDataPtr;
117 OpError(boost::shared_ptr<CommonData> &common_data_ptr)
119 std::fill(doEntities.begin(), doEntities.end(),
false);
120 doEntities[MBVERTEX] =
true;
125 const int nb_integration_pts = getGaussPts().size2();
126 auto t_w = getFTensor0IntegrationWeight();
128 commonDataPtr->approxVals));
129 auto t_grad_val = getFTensor1FromMat<SPACE_DIM>(
132 auto t_hessian_val = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
133 *(commonDataPtr)->approxHessianVals);
136 auto t_inv_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
137 *(commonDataPtr->invJacPtr));
138 auto t_coords = getFTensor1CoordsAtGaussPts();
147 const double volume = getMeasure();
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};
216 CHKERR mField.getInterface(simpleInterface);
217 CHKERR simpleInterface->getOptions();
218 CHKERR simpleInterface->loadFile();
230 constexpr
int order = 4;
232 CHKERR simpleInterface->setUp();
242 auto rule = [](
int,
int,
int p) ->
int {
return 2 * p; };
263 MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
266 CHKERR KSPSetFromOptions(solver);
269 auto dm = simpleInterface->getDM();
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>();
295 mField.get_comm(), (!mField.get_comm_rank()) ? 3 : 0, 3);
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);
397 int main(
int argc,
char *argv[]) {
405 DMType dm_name =
"DMMOFEM";