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:
boost::shared_ptr<MatrixDouble> invJacPtr;
boost::shared_ptr<VectorDouble> approxVals;
boost::shared_ptr<MatrixDouble> approxGradVals;
boost::shared_ptr<MatrixDouble> approxHessianVals;
};
};
boost::shared_ptr<CommonData> commonDataPtr;
OpError(boost::shared_ptr<CommonData> &common_data_ptr)
std::fill(doEntities.begin(), doEntities.end(), false);
doEntities[MBVERTEX] = true;
}
const int nb_integration_pts = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
commonDataPtr->approxVals));
auto t_grad_val = getFTensor1FromMat<SPACE_DIM>(
*(commonDataPtr
->approxGradVals));
auto t_hessian_val = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
*(commonDataPtr)->approxHessianVals);
auto t_inv_jac = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
*(commonDataPtr->invJacPtr));
auto t_coords = getFTensor1CoordsAtGaussPts();
const double volume = getMeasure();
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);
}
};
}
CHKERR mField.getInterface(simpleInterface);
CHKERR simpleInterface->getOptions();
CHKERR simpleInterface->loadFile();
}
CHKERR simpleInterface->setUp();
}
auto rule = [](
int,
int,
int p) ->
int {
return 2 * p; };
}
MOFEM_LOG(
"WORLD", Sev::inform) <<
"Solve problem";
CHKERR KSPSetFromOptions(solver);
auto dm = simpleInterface->getDM();
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>();
mField.get_comm(), (!mField.get_comm_rank()) ? 3 : 0, 3);
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";
}
}