v0.10.0
Functions | Variables
testing_jacobian_of_hook_element.cpp File Reference
#include <BasicFiniteElements.hpp>

Go to the source code of this file.

Functions

int main (int argc, char *argv[])
 

Variables

static char help [] = "\n"
 

Function Documentation

◆ main()

int main ( int  argc,
char *  argv[] 
)
Examples
testing_jacobian_of_hook_element.cpp.

Definition at line 29 of file testing_jacobian_of_hook_element.cpp.

29  {
30 
31  // Initialize MoFEM
32  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
33 
34  // Create mesh database
35  moab::Core mb_instance; // create database
36  moab::Interface &moab = mb_instance; // create interface to database
37 
38  try {
39  // Create MoFEM database and link it to MoAB
40  MoFEM::Core core(moab);
41  MoFEM::Interface &m_field = core;
42 
43  PetscBool ale = PETSC_FALSE;
44  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-ale", &ale, PETSC_NULL);
45  PetscBool test_jacobian = PETSC_FALSE;
46  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_jacobian", &test_jacobian,
47  PETSC_NULL);
48 
49  CHKERR DMRegister_MoFEM("DMMOFEM");
50 
51  Simple *si = m_field.getInterface<MoFEM::Simple>();
52 
53  CHKERR si->getOptions();
54  CHKERR si->loadFile();
56  const int order = 2;
57  CHKERR si->setFieldOrder("x", order);
58 
59  if (ale == PETSC_TRUE) {
61  CHKERR si->setFieldOrder("X", 2);
62  }
63 
64  CHKERR si->setUp();
65 
66  // create DM
67  DM dm;
68  CHKERR si->getDM(&dm);
69 
70  // Projection on "x" field
71  {
72  Projection10NodeCoordsOnField ent_method(m_field, "x");
73  CHKERR m_field.loop_dofs("x", ent_method);
74  // CHKERR m_field.getInterface<FieldBlas>()->fieldScale(2, "x");
75  }
76 
77  // Project coordinates on "X" field
78  if (ale == PETSC_TRUE) {
79  Projection10NodeCoordsOnField ent_method(m_field, "X");
80  CHKERR m_field.loop_dofs("X", ent_method);
81  // CHKERR m_field.getInterface<FieldBlas>()->fieldScale(2, "X");
82  }
83 
84  boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
86  boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
88  struct VolRule {
89  int operator()(int, int, int) const { return 2 * (order - 1); }
90  };
91  fe_lhs_ptr->getRuleHook = VolRule();
92  fe_rhs_ptr->getRuleHook = VolRule();
93 
94  CHKERR DMMoFEMSNESSetJacobian(dm, si->getDomainFEName(), fe_lhs_ptr,
95  nullptr, nullptr);
96  CHKERR DMMoFEMSNESSetFunction(dm, si->getDomainFEName(), fe_rhs_ptr,
97  nullptr, nullptr);
98 
100  boost::shared_ptr<map<int, BlockData>> block_sets_ptr =
101  boost::make_shared<map<int, BlockData>>();
102  (*block_sets_ptr)[0].iD = 0;
103  (*block_sets_ptr)[0].E = 1;
104  (*block_sets_ptr)[0].PoissonRatio = 0.25;
106  si->getDomainFEName(), 3, (*block_sets_ptr)[0].tEts);
107 
108  CHKERR HookeElement::setOperators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
109  "x", "X", ale, false);
110 
111  Vec x, f;
112  CHKERR DMCreateGlobalVector(dm, &x);
113  CHKERR VecDuplicate(x, &f);
114  CHKERR DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_FORWARD);
115 
116  // CHKERR VecDuplicate(x, &dx);
117  // PetscRandom rctx;
118  // PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
119  // VecSetRandom(dx, rctx);
120  // PetscRandomDestroy(&rctx);
121  // CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);
122 
123  Mat A, fdA;
124  CHKERR DMCreateMatrix(dm, &A);
125  CHKERR MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &fdA);
126 
127  if (test_jacobian == PETSC_TRUE) {
128  char testing_options[] =
129  "-snes_test_jacobian -snes_test_jacobian_display "
130  "-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 -snes_max_it 1 "
131  "-pc_type none";
132  CHKERR PetscOptionsInsertString(NULL, testing_options);
133  } else {
134  char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
135  "-snes_rtol 0 -snes_max_it 1 -pc_type none";
136  CHKERR PetscOptionsInsertString(NULL, testing_options);
137  }
138 
139  SNES snes;
140  CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
141  MoFEM::SnesCtx *snes_ctx;
142  CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
143  CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
144  CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
145  CHKERR SNESSetFromOptions(snes);
146 
147  CHKERR SNESSolve(snes, NULL, x);
148 
149  if (test_jacobian == PETSC_FALSE) {
150  double nrm_A0;
151  CHKERR MatNorm(A, NORM_INFINITY, &nrm_A0);
152 
153  char testing_options_fd[] = "-snes_fd";
154  CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
155 
156  CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
157  CHKERR SNESSetJacobian(snes, fdA, fdA, SnesMat, snes_ctx);
158  CHKERR SNESSetFromOptions(snes);
159 
160  CHKERR SNESSolve(snes, NULL, x);
161  CHKERR MatAXPY(A, -1, fdA, SUBSET_NONZERO_PATTERN);
162 
163  double nrm_A;
164  CHKERR MatNorm(A, NORM_INFINITY, &nrm_A);
165  PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
166  nrm_A / nrm_A0);
167  nrm_A /= nrm_A0;
168 
169  const double tol = 1e-5;
170  if (nrm_A > tol) {
171  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
172  "Difference between hand-calculated tangent matrix and finite "
173  "difference matrix is too big");
174  }
175  }
176 
177  CHKERR VecDestroy(&x);
178  CHKERR VecDestroy(&f);
179  CHKERR MatDestroy(&A);
180  CHKERR MatDestroy(&fdA);
181  CHKERR SNESDestroy(&snes);
182 
183  // destroy DM
184  CHKERR DMDestroy(&dm);
185  }
186  CATCH_ERRORS;
187 
188  // finish work cleaning memory, getting statistics, etc
190 
191  return 0;
192 }
Set integration rule.
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:672
PetscErrorCode DMMoFEMSNESSetFunction(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES residual evaluation function
Definition: DMMMoFEM.cpp:637
Deprecated interface functions.
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
static MoFEMErrorCode setOperators(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, const EntityType type=MBTET, boost::shared_ptr< DataAtIntegrationPts > data_at_pts=nullptr)
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:706
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
CoreTmp< 0 > Core
Definition: Core.hpp:1129
int operator()(int, int, int) const
const FTensor::Tensor2< T, Dim, Dim > Vec
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:445
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.
Definition: SnesCtx.cpp:126
static char help[]
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
Projection of edge entities with one mid-node on hierarchical basis.
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
constexpr int order
double tol
NonlinearElasticElement::BlockData BlockData
Definition: elasticity.cpp:89
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:212
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMMoFEM.cpp:953
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:27
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
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.
Definition: Simple.cpp:269
PetscErrorCode DMMoFEMSNESSetJacobian(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set SNES Jacobian evaluation function
Definition: DMMMoFEM.cpp:678
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.
#define CHKERR
Inline error check.
Definition: definitions.h:604
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:281
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:198
Core (interface) class.
Definition: Core.hpp:77
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
continuous field
Definition: definitions.h:177
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:482
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.
Definition: SnesCtx.cpp:17
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100

Variable Documentation

◆ help

char help[] = "\n"
static