Testing implementation of Hook element by verifying tangent stiffness matrix. Test like this is an example of how to verify the implementation of Jacobian.
#include <BasicFiniteElements.hpp>
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) {
CHKERR si->setFieldOrder(
"X", 2);
}
DM dm;
{
Projection10NodeCoordsOnField ent_method(m_field, "x");
}
if (ale == PETSC_TRUE) {
Projection10NodeCoordsOnField ent_method(m_field, "X");
}
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;
si->getDomainFEName(), 3, (*block_sets_ptr)[0].tEts);
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(
new OpCalculateVectorFieldValues<3>(x_field, mat_coords_ptr));
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(
new OpCalculateVectorFieldValues<3>(X_field, mat_coords_ptr));
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(
new OpCalculateVectorFieldValues<3>(x_field, mat_coords_ptr));
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(
new OpCalculateVectorFieldValues<3>(X_field, mat_coords_ptr));
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;
}
EntitiesFieldData::EntData EntData
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_ATOM_TEST_INVALID
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
NonlinearElasticElement::BlockData BlockData
VolumeElementForcesAndSourcesCore::UserDataOperator VolUserDataOperator
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
virtual MoFEMErrorCode get_finite_element_entities_by_dimension(const std::string name, int dim, Range &ents) const =0
get entities in the finite element by dimension
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
FTensor::Index< 'i', SPACE_DIM > i
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
UBlasMatrix< double > MatrixDouble
implementation of Data Operators for Forces and Sources
PetscErrorCode SnesMat(SNES snes, Vec x, Mat A, Mat B, void *ctx)
This is MoFEM implementation for the left hand side (tangent matrix) evaluation in SNES solver.
PetscErrorCode SnesRhs(SNES snes, Vec x, Vec f, void *ctx)
This is MoFEM implementation for the right hand side (residual vector) evaluation in SNES solver.
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
DeprecatedCoreInterface Interface
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.
Deprecated interface functions.
Simple interface for fast problem set-up.
Interface for nonlinear (SNES) solver.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
int main(int argc, char *argv[])