Testing implementation of Hook element by verifying tangent stiffness matrix. Test like this is an example of how to verify the implementation of Jacobian.
using namespace boost::numeric;
static char help[] =
"\n";
template <bool ALE>
boost::shared_ptr<MatrixDouble> matCoordsPtr;
boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
boost::shared_ptr<MatrixDouble> mat_coords_ptr,
boost::shared_ptr<VectorDouble> density_at_pts,
boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr)
matCoordsPtr(mat_coords_ptr), rhoAtGaussPtsPtr(density_at_pts),
rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts_ptr) {}
if (row_type != MBVERTEX)
const int nb_integration_pts = getGaussPts().size2();
rhoAtGaussPtsPtr->resize(nb_integration_pts, false);
rhoAtGaussPtsPtr->clear();
rhoGradAtGaussPtsPtr->resize(3, nb_integration_pts, false);
rhoGradAtGaussPtsPtr->clear();
auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
if (ALE)
coords = *matCoordsPtr;
else
coords = trans(getCoordsAtGaussPts());
auto t_coords = getFTensor1FromMat<3>(coords);
for (int gg = 0; gg != nb_integration_pts; ++gg) {
t_rho = 1 + t_coords(
i) * t_coords(
i);
t_grad_rho(
i) = 2 * t_coords(
i);
++t_rho;
++t_coords;
++t_grad_rho;
}
}
};
int main(
int argc,
char *argv[]) {
try {
PetscBool ale = PETSC_FALSE;
PetscBool test_jacobian = PETSC_FALSE;
PETSC_NULL);
if (ale == PETSC_TRUE) {
}
DM dm;
{
}
if (ale == PETSC_TRUE) {
}
boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
int operator()(
int,
int,
int)
const {
return 2 * (
order - 1); }
};
fe_lhs_ptr->getRuleHook =
VolRule();
fe_rhs_ptr->getRuleHook =
VolRule();
nullptr, nullptr);
nullptr, nullptr);
boost::shared_ptr<map<int, BlockData>> block_sets_ptr =
boost::make_shared<map<int, BlockData>>();
(*block_sets_ptr)[0].
iD = 0;
(*block_sets_ptr)[0].E = 1;
(*block_sets_ptr)[0].PoissonRatio = 0.25;
const double rho_n = 2.0;
const double rho_0 = 0.5;
auto my_operators = [&](boost::shared_ptr<ForcesAndSourcesCore> &fe_lhs_ptr,
boost::shared_ptr<ForcesAndSourcesCore> &fe_rhs_ptr,
boost::shared_ptr<map<int, BlockData>>
&block_sets_ptr,
const std::string x_field,
const std::string X_field, const bool ale,
const bool field_disp) {
boost::shared_ptr<HookeElement::DataAtIntegrationPts> data_at_pts(
new HookeElement::DataAtIntegrationPts());
boost::shared_ptr<MatrixDouble> mat_coords_ptr =
boost::make_shared<MatrixDouble>();
boost::shared_ptr<VectorDouble> rho_at_gauss_pts_ptr =
boost::make_shared<VectorDouble>();
boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr =
boost::make_shared<MatrixDouble>();
if (fe_lhs_ptr) {
if (ale == PETSC_FALSE) {
fe_lhs_ptr->getOpPtrVector().push_back(
x_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
rho_grad_at_gauss_pts_ptr));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStiffnessScaledByDensityField(
x_field, x_field, block_sets_ptr, data_at_pts,
rho_at_gauss_pts_ptr, rho_n, rho_0));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpLhs_dx_dx<1>(x_field, x_field, data_at_pts));
} else {
fe_lhs_ptr->getOpPtrVector().push_back(
data_at_pts->HMat));
fe_lhs_ptr->getOpPtrVector().push_back(
X_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
rho_grad_at_gauss_pts_ptr));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStiffnessScaledByDensityField(
x_field, x_field, block_sets_ptr, data_at_pts,
rho_at_gauss_pts_ptr, rho_n, rho_0));
fe_lhs_ptr->getOpPtrVector().push_back(
data_at_pts->hMat));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStrainAle(x_field, x_field,
data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStress<1>(x_field, x_field,
data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dx_dx<1>(x_field, x_field,
data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dx_dX<1>(x_field, X_field,
data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateEnergy(X_field, X_field,
data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dX_dX<1>(X_field, X_field,
data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAleLhsPre_dX_dx<1>(X_field, x_field,
data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAleLhs_dX_dx(X_field, x_field, data_at_pts));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAleLhsWithDensity_dx_dX(
x_field, X_field, data_at_pts, rho_at_gauss_pts_ptr,
rho_grad_at_gauss_pts_ptr, rho_n, rho_0));
fe_lhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAleLhsWithDensity_dX_dX(
X_field, X_field, data_at_pts, rho_at_gauss_pts_ptr,
rho_grad_at_gauss_pts_ptr, rho_n, rho_0));
}
}
if (fe_rhs_ptr) {
if (ale == PETSC_FALSE) {
fe_rhs_ptr->getOpPtrVector().push_back(
data_at_pts->hMat));
fe_rhs_ptr->getOpPtrVector().push_back(
x_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
rho_grad_at_gauss_pts_ptr));
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStiffnessScaledByDensityField(
x_field, x_field, block_sets_ptr, data_at_pts,
rho_at_gauss_pts_ptr, rho_n, rho_0));
if (field_disp) {
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStrain<1>(x_field, x_field,
data_at_pts));
} else {
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStrain<0>(x_field, x_field,
data_at_pts));
}
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStress<1>(x_field, x_field,
data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpRhs_dx(x_field, x_field, data_at_pts));
} else {
fe_rhs_ptr->getOpPtrVector().push_back(
data_at_pts->HMat));
fe_rhs_ptr->getOpPtrVector().push_back(
X_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
rho_grad_at_gauss_pts_ptr));
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStiffnessScaledByDensityField(
x_field, x_field, block_sets_ptr, data_at_pts,
rho_at_gauss_pts_ptr, rho_n, rho_0));
fe_rhs_ptr->getOpPtrVector().push_back(
data_at_pts->hMat));
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStrainAle(x_field, x_field,
data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateStress<1>(x_field, x_field,
data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAleRhs_dx(x_field, x_field, data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateEnergy(X_field, X_field,
data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
data_at_pts));
fe_rhs_ptr->getOpPtrVector().push_back(
new HookeElement::OpAleRhs_dX(X_field, X_field, data_at_pts));
}
}
};
CHKERR my_operators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
"x",
"X", ale,
false);
CHKERR DMCreateGlobalVector(dm, &x);
CHKERR MatDuplicate(
A, MAT_DO_NOT_COPY_VALUES, &fdA);
if (test_jacobian == PETSC_TRUE) {
char testing_options[] =
"-snes_test_jacobian -snes_test_jacobian_display "
"-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 -snes_max_it 1 "
"-pc_type none";
CHKERR PetscOptionsInsertString(NULL, testing_options);
} else {
char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
"-snes_rtol 0 -snes_max_it 1 -pc_type none";
CHKERR PetscOptionsInsertString(NULL, testing_options);
}
SNES snes;
CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
CHKERR SNESSetFromOptions(snes);
CHKERR SNESSolve(snes, NULL, x);
if (test_jacobian == PETSC_FALSE) {
double nrm_A0;
CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A0);
char testing_options_fd[] = "-snes_fd";
CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
CHKERR SNESSetFromOptions(snes);
CHKERR SNESSolve(snes, NULL, x);
CHKERR MatAXPY(
A, -1, fdA, SUBSET_NONZERO_PATTERN);
double nrm_A;
CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A);
PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
nrm_A / nrm_A0);
nrm_A /= nrm_A0;
"Difference between hand-calculated tangent matrix and finite "
"difference matrix is too big");
}
}
}
return 0;
}