v0.14.0
Loading...
Searching...
No Matches
testing_jacobian_of_hook_scaled_with_density_element.cpp
Go to the documentation of this file.
1/** \file testing_jacobian_of_hook_scaled_with_density_element.cpp
2 * \example testing_jacobian_of_hook_scaled_with_density_element.cpp
3
4Testing implementation of Hook element by verifying tangent stiffness matrix.
5Test like this is an example of how to verify the implementation of Jacobian.
6
7*/
8
9
10
12
13using namespace boost::numeric;
14using namespace MoFEM;
15
16static char help[] = "\n";
17
18template <bool ALE>
19struct OpGetDensityField : public HookeElement::VolUserDataOperator {
20
21 boost::shared_ptr<MatrixDouble> matCoordsPtr;
22 boost::shared_ptr<VectorDouble> rhoAtGaussPtsPtr;
23 boost::shared_ptr<MatrixDouble> rhoGradAtGaussPtsPtr;
24 OpGetDensityField(const std::string row_field,
25 boost::shared_ptr<MatrixDouble> mat_coords_ptr,
26 boost::shared_ptr<VectorDouble> density_at_pts,
27 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr)
28 : HookeElement::VolUserDataOperator(row_field, OPROW),
29 matCoordsPtr(mat_coords_ptr), rhoAtGaussPtsPtr(density_at_pts),
30 rhoGradAtGaussPtsPtr(rho_grad_at_gauss_pts_ptr) {}
31
32 MoFEMErrorCode doWork(int row_side, EntityType row_type,
33 HookeElement::EntData &row_data) {
35 if (row_type != MBVERTEX)
37 // get number of integration points
38 const int nb_integration_pts = getGaussPts().size2();
39 rhoAtGaussPtsPtr->resize(nb_integration_pts, false);
40 rhoAtGaussPtsPtr->clear();
41 rhoGradAtGaussPtsPtr->resize(3, nb_integration_pts, false);
42 rhoGradAtGaussPtsPtr->clear();
43
45 auto t_grad_rho = getFTensor1FromMat<3>(*rhoGradAtGaussPtsPtr);
46
47 MatrixDouble coords;
48 if (ALE)
49 coords = *matCoordsPtr;
50 else
51 coords = trans(getCoordsAtGaussPts()); // because the size is (nb_gg,3)
52
54 auto t_coords = getFTensor1FromMat<3>(coords);
55
56 for (int gg = 0; gg != nb_integration_pts; ++gg) {
57 t_rho = 1 + t_coords(i) * t_coords(i); // density is equal to
58 // 1+x^2+y^2+z^2 (x,y,z coordines)
59 t_grad_rho(i) = 2 * t_coords(i);
60
61 ++t_rho;
62 ++t_coords;
63 ++t_grad_rho;
64 }
65
67 }
68};
69
70int main(int argc, char *argv[]) {
71
72 // Initialize MoFEM
73 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
74
75 // Create mesh database
76 moab::Core mb_instance; // create database
77 moab::Interface &moab = mb_instance; // create interface to database
78
79 try {
80 // Create MoFEM database and link it to MoAB
81 MoFEM::Core core(moab);
82 MoFEM::Interface &m_field = core;
83
84 PetscBool ale = PETSC_FALSE;
85 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-ale", &ale, PETSC_NULL);
86 PetscBool test_jacobian = PETSC_FALSE;
87 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test_jacobian", &test_jacobian,
88 PETSC_NULL);
89
90 CHKERR DMRegister_MoFEM("DMMOFEM");
91
92 Simple *si = m_field.getInterface<MoFEM::Simple>();
93
94 CHKERR si->getOptions();
95 CHKERR si->loadFile();
97 const int order = 2;
98 CHKERR si->setFieldOrder("x", order);
99
100 if (ale == PETSC_TRUE) {
102 CHKERR si->setFieldOrder("X", 2);
103 }
104 CHKERR si->setUp();
105
106 // create DM
107 DM dm;
108 CHKERR si->getDM(&dm);
109
110 // Projection on "x" field
111 {
112 Projection10NodeCoordsOnField ent_method(m_field, "x");
113 CHKERR m_field.loop_dofs("x", ent_method);
114 // CHKERR m_field.getInterface<FieldBlas>()->fieldScale(1.5, "x");
115 }
116
117 // Project coordinates on "X" field
118 if (ale == PETSC_TRUE) {
119 Projection10NodeCoordsOnField ent_method(m_field, "X");
120 CHKERR m_field.loop_dofs("X", ent_method);
121 // CHKERR m_field.getInterface<FieldBlas>()->fieldScale(1.5, "X");
122 }
123
124 boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
126 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
128 struct VolRule {
129 int operator()(int, int, int) const { return 2 * (order - 1); }
130 };
131 fe_lhs_ptr->getRuleHook = VolRule();
132 fe_rhs_ptr->getRuleHook = VolRule();
133
134 CHKERR DMMoFEMSNESSetJacobian(dm, si->getDomainFEName(), fe_lhs_ptr,
135 nullptr, nullptr);
136 CHKERR DMMoFEMSNESSetFunction(dm, si->getDomainFEName(), fe_rhs_ptr,
137 nullptr, nullptr);
138
140 boost::shared_ptr<map<int, BlockData>> block_sets_ptr =
141 boost::make_shared<map<int, BlockData>>();
142 (*block_sets_ptr)[0].iD = 0;
143 (*block_sets_ptr)[0].E = 1;
144 (*block_sets_ptr)[0].PoissonRatio = 0.25;
146 si->getDomainFEName(), 3, (*block_sets_ptr)[0].tEts);
147
148 // Parameters for density
149 const double rho_n = 2.0;
150 const double rho_0 = 0.5;
151
152 auto my_operators = [&](boost::shared_ptr<ForcesAndSourcesCore> &fe_lhs_ptr,
153 boost::shared_ptr<ForcesAndSourcesCore> &fe_rhs_ptr,
154 boost::shared_ptr<map<int, BlockData>>
155 &block_sets_ptr,
156 const std::string x_field,
157 const std::string X_field, const bool ale,
158 const bool field_disp) {
160
161 boost::shared_ptr<HookeElement::DataAtIntegrationPts> data_at_pts(
162 new HookeElement::DataAtIntegrationPts());
163 boost::shared_ptr<MatrixDouble> mat_coords_ptr =
164 boost::make_shared<MatrixDouble>();
165 boost::shared_ptr<VectorDouble> rho_at_gauss_pts_ptr =
166 boost::make_shared<VectorDouble>();
167 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr =
168 boost::make_shared<MatrixDouble>();
169
170 if (fe_lhs_ptr) {
171 if (ale == PETSC_FALSE) {
172 fe_lhs_ptr->getOpPtrVector().push_back(
173 new OpCalculateVectorFieldValues<3>(x_field, mat_coords_ptr));
174 fe_lhs_ptr->getOpPtrVector().push_back(new OpGetDensityField<false>(
175 x_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
176 rho_grad_at_gauss_pts_ptr));
177 fe_lhs_ptr->getOpPtrVector().push_back(
178 new HookeElement::OpCalculateStiffnessScaledByDensityField(
179 x_field, x_field, block_sets_ptr, data_at_pts,
180 rho_at_gauss_pts_ptr, rho_n, rho_0));
181 fe_lhs_ptr->getOpPtrVector().push_back(
182 new HookeElement::OpLhs_dx_dx<1>(x_field, x_field, data_at_pts));
183 } else {
184 fe_lhs_ptr->getOpPtrVector().push_back(
186 data_at_pts->HMat));
187 fe_lhs_ptr->getOpPtrVector().push_back(
188 new OpCalculateVectorFieldValues<3>(X_field, mat_coords_ptr));
189 fe_lhs_ptr->getOpPtrVector().push_back(new OpGetDensityField<true>(
190 X_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
191 rho_grad_at_gauss_pts_ptr));
192 fe_lhs_ptr->getOpPtrVector().push_back(
193 new HookeElement::OpCalculateStiffnessScaledByDensityField(
194 x_field, x_field, block_sets_ptr, data_at_pts,
195 rho_at_gauss_pts_ptr, rho_n, rho_0));
196 fe_lhs_ptr->getOpPtrVector().push_back(
198 data_at_pts->hMat));
199 fe_lhs_ptr->getOpPtrVector().push_back(
200 new HookeElement::OpCalculateStrainAle(x_field, x_field,
201 data_at_pts));
202 fe_lhs_ptr->getOpPtrVector().push_back(
203 new HookeElement::OpCalculateStress<1>(x_field, x_field,
204 data_at_pts)); // FIXME:
205 fe_lhs_ptr->getOpPtrVector().push_back(
206 new HookeElement::OpAleLhs_dx_dx<1>(x_field, x_field,
207 data_at_pts));
208 fe_lhs_ptr->getOpPtrVector().push_back(
209 new HookeElement::OpAleLhs_dx_dX<1>(x_field, X_field,
210 data_at_pts));
211 fe_lhs_ptr->getOpPtrVector().push_back(
212 new HookeElement::OpCalculateEnergy(X_field, X_field,
213 data_at_pts));
214 fe_lhs_ptr->getOpPtrVector().push_back(
215 new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
216 data_at_pts));
217 fe_lhs_ptr->getOpPtrVector().push_back(
218 new HookeElement::OpAleLhs_dX_dX<1>(X_field, X_field,
219 data_at_pts));
220 fe_lhs_ptr->getOpPtrVector().push_back(
221 new HookeElement::OpAleLhsPre_dX_dx<1>(X_field, x_field,
222 data_at_pts));
223 fe_lhs_ptr->getOpPtrVector().push_back(
224 new HookeElement::OpAleLhs_dX_dx(X_field, x_field, data_at_pts));
225 fe_lhs_ptr->getOpPtrVector().push_back(
226 new HookeElement::OpAleLhsWithDensity_dx_dX(
227 x_field, X_field, data_at_pts, rho_at_gauss_pts_ptr,
228 rho_grad_at_gauss_pts_ptr, rho_n, rho_0));
229 fe_lhs_ptr->getOpPtrVector().push_back(
230 new HookeElement::OpAleLhsWithDensity_dX_dX(
231 X_field, X_field, data_at_pts, rho_at_gauss_pts_ptr,
232 rho_grad_at_gauss_pts_ptr, rho_n, rho_0));
233 }
234 }
235
236 if (fe_rhs_ptr) {
237
238 if (ale == PETSC_FALSE) {
239 fe_rhs_ptr->getOpPtrVector().push_back(
241 data_at_pts->hMat));
242 fe_rhs_ptr->getOpPtrVector().push_back(
243 new OpCalculateVectorFieldValues<3>(x_field, mat_coords_ptr));
244 fe_rhs_ptr->getOpPtrVector().push_back(new OpGetDensityField<false>(
245 x_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
246 rho_grad_at_gauss_pts_ptr));
247 fe_rhs_ptr->getOpPtrVector().push_back(
248 new HookeElement::OpCalculateStiffnessScaledByDensityField(
249 x_field, x_field, block_sets_ptr, data_at_pts,
250 rho_at_gauss_pts_ptr, rho_n, rho_0));
251 if (field_disp) {
252 fe_rhs_ptr->getOpPtrVector().push_back(
253 new HookeElement::OpCalculateStrain<1>(x_field, x_field,
254 data_at_pts));
255 } else {
256 fe_rhs_ptr->getOpPtrVector().push_back(
257 new HookeElement::OpCalculateStrain<0>(x_field, x_field,
258 data_at_pts));
259 }
260 fe_rhs_ptr->getOpPtrVector().push_back(
261 new HookeElement::OpCalculateStress<1>(x_field, x_field,
262 data_at_pts));
263 fe_rhs_ptr->getOpPtrVector().push_back(
264 new HookeElement::OpRhs_dx(x_field, x_field, data_at_pts));
265 } else {
266 fe_rhs_ptr->getOpPtrVector().push_back(
268 data_at_pts->HMat));
269 fe_rhs_ptr->getOpPtrVector().push_back(
270 new OpCalculateVectorFieldValues<3>(X_field, mat_coords_ptr));
271 fe_rhs_ptr->getOpPtrVector().push_back(new OpGetDensityField<true>(
272 X_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
273 rho_grad_at_gauss_pts_ptr));
274 fe_rhs_ptr->getOpPtrVector().push_back(
275 new HookeElement::OpCalculateStiffnessScaledByDensityField(
276 x_field, x_field, block_sets_ptr, data_at_pts,
277 rho_at_gauss_pts_ptr, rho_n, rho_0));
278 fe_rhs_ptr->getOpPtrVector().push_back(
280 data_at_pts->hMat));
281 fe_rhs_ptr->getOpPtrVector().push_back(
282 new HookeElement::OpCalculateStrainAle(x_field, x_field,
283 data_at_pts));
284 fe_rhs_ptr->getOpPtrVector().push_back(
285 new HookeElement::OpCalculateStress<1>(x_field, x_field,
286 data_at_pts));
287 fe_rhs_ptr->getOpPtrVector().push_back(
288 new HookeElement::OpAleRhs_dx(x_field, x_field, data_at_pts));
289 fe_rhs_ptr->getOpPtrVector().push_back(
290 new HookeElement::OpCalculateEnergy(X_field, X_field,
291 data_at_pts));
292 fe_rhs_ptr->getOpPtrVector().push_back(
293 new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
294 data_at_pts));
295 fe_rhs_ptr->getOpPtrVector().push_back(
296 new HookeElement::OpAleRhs_dX(X_field, X_field, data_at_pts));
297 }
298 }
300 };
301 CHKERR my_operators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr, "x", "X", ale,
302 false);
303 Vec x, f;
304 CHKERR DMCreateGlobalVector(dm, &x);
305 CHKERR VecDuplicate(x, &f);
306 CHKERR DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_FORWARD);
307
308 // CHKERR VecDuplicate(x, &dx);
309 // PetscRandom rctx;
310 // PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
311 // VecSetRandom(dx, rctx);
312 // PetscRandomDestroy(&rctx);
313 // CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);
314
315 Mat A, fdA;
316 CHKERR DMCreateMatrix(dm, &A);
317 CHKERR MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &fdA);
318
319 if (test_jacobian == PETSC_TRUE) {
320 char testing_options[] =
321 "-snes_test_jacobian -snes_test_jacobian_display "
322 "-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 -snes_max_it 1 "
323 "-pc_type none";
324 CHKERR PetscOptionsInsertString(NULL, testing_options);
325 } else {
326 char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
327 "-snes_rtol 0 -snes_max_it 1 -pc_type none";
328 CHKERR PetscOptionsInsertString(NULL, testing_options);
329 }
330
331 SNES snes;
332 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
333 MoFEM::SnesCtx *snes_ctx;
334 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
335 CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
336 CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
337 CHKERR SNESSetFromOptions(snes);
338
339 CHKERR SNESSolve(snes, NULL, x);
340
341 if (test_jacobian == PETSC_FALSE) {
342 double nrm_A0;
343 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A0);
344
345 char testing_options_fd[] = "-snes_fd";
346 CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
347
348 CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
349 CHKERR SNESSetJacobian(snes, fdA, fdA, SnesMat, snes_ctx);
350 CHKERR SNESSetFromOptions(snes);
351
352 CHKERR SNESSolve(snes, NULL, x);
353 CHKERR MatAXPY(A, -1, fdA, SUBSET_NONZERO_PATTERN);
354
355 double nrm_A;
356 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A);
357 PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
358 nrm_A / nrm_A0);
359 nrm_A /= nrm_A0;
360
361 const double tol = 1e-5;
362 if (nrm_A > tol) {
363 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
364 "Difference between hand-calculated tangent matrix and finite "
365 "difference matrix is too big");
366 }
367 }
368
369 CHKERR VecDestroy(&x);
370 CHKERR VecDestroy(&f);
371 CHKERR MatDestroy(&A);
372 CHKERR MatDestroy(&fdA);
373 CHKERR SNESDestroy(&snes);
374
375 // destroy DM
376 CHKERR DMDestroy(&dm);
377 }
379
380 // finish work cleaning memory, getting statistics, etc
382
383 return 0;
384}
int main()
Definition: adol-c_atom.cpp:46
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
constexpr int order
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: DMMoFEM.cpp:704
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
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: DMMoFEM.cpp:745
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1080
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
double tol
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
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:136
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:27
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.
Definition: Templates.hpp:135
constexpr AssemblyType A
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
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.
Definition: Simple.hpp:27
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
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:264
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:341
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
data for calculation heat conductivity and heat capacity elements
MoFEMErrorCode doWork(int row_side, EntityType row_type, HookeElement::EntData &row_data)
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
OpGetDensityField(const std::string row_field, boost::shared_ptr< MatrixDouble > mat_coords_ptr, boost::shared_ptr< VectorDouble > density_at_pts, boost::shared_ptr< MatrixDouble > rho_grad_at_gauss_pts_ptr)
Set integration rule.
int operator()(int, int, int) const