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> mat_coords_ptr,
boost::shared_ptr<VectorDouble> density_at_pts,
boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr)
: HookeElement::VolUserDataOperator(row_field, OPROW),
HookeElement::EntData &row_data) {
if (row_type != MBVERTEX)
const int nb_integration_pts = getGaussPts().size2();
if (ALE)
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[]) {
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
try {
PetscBool ale = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL,
"",
"-ale", &ale, PETSC_NULL);
PetscBool test_jacobian = PETSC_FALSE;
CHKERR PetscOptionsGetBool(PETSC_NULL,
"",
"-test_jacobian", &test_jacobian,
PETSC_NULL);
CHKERR DMRegister_MoFEM(
"DMMOFEM");
if (ale == PETSC_TRUE) {
}
DM dm;
{
}
if (ale == PETSC_TRUE) {
}
boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
};
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 SNESSetFunction(snes, f, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes,
A,
A, SnesMat, snes_ctx);
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 SNESSetFunction(snes, f, SnesRhs, snes_ctx);
CHKERR SNESSetJacobian(snes, fdA, fdA, SnesMat, snes_ctx);
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;
}
#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 ...
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 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
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
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.
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Projection of edge entities with one mid-node on hierarchical basis.
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.
const std::string getDomainFEName() const
Get the Domain FE Name.
Interface for nonlinear (SNES) solver.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
data for calculation heat conductivity and heat capacity elements
MoFEMErrorCode doWork(int row_side, EntityType row_type, HookeElement::EntData &row_data)
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > matCoordsPtr
int operator()(int, int, int) const