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";
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;
CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
"x", "X", ale, false);
Vec x, f;
CHKERR DMCreateGlobalVector(dm, &x);
CHKERR DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_FORWARD);
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 DMMoFEMGetSnesCtx(dm, &snes_ctx);
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 .
@ MOFEM_ATOM_TEST_INVALID
#define CHKERR
Inline error check.
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.
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.
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
int operator()(int, int, int) const