v0.16.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 mofem/users_modules/basic_finite_elements/atom_tests/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 =
46 getFTensor1FromMat<3, -1, CoeffsByGauss>(*rhoGradAtGaussPtsPtr);
47
48 FTensor::Index<'i', 3> i;
49 auto set_density = [&](auto t_coords) {
50 for (int gg = 0; gg != nb_integration_pts; ++gg) {
51 t_rho = 1 + t_coords(i) * t_coords(i); // 1+x^2+y^2+z^2
52 t_grad_rho(i) = 2 * t_coords(i);
53
54 ++t_rho;
55 ++t_coords;
56 ++t_grad_rho;
57 }
58 };
59
60 if constexpr (ALE) {
61 set_density(getFTensor1FromMat<3>(*matCoordsPtr));
62 } else {
63 MatrixDouble coords =
64 trans(getCoordsAtGaussPts()); // because the size is (nb_gg,3)
65 set_density(getFTensor1FromMat<3, -1, CoeffsByGauss>(coords));
66 }
67
69 }
70};
71
72int main(int argc, char *argv[]) {
73
74 // Initialize MoFEM
75 MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
76
77 // Create mesh database
78 moab::Core mb_instance; // create database
79 moab::Interface &moab = mb_instance; // create interface to database
80
81 try {
82 // Create MoFEM database and link it to MoAB
83 MoFEM::Core core(moab);
84 MoFEM::Interface &m_field = core;
85
86 PetscBool ale = PETSC_FALSE;
87 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-ale", &ale, PETSC_NULLPTR);
88 PetscBool test_jacobian = PETSC_FALSE;
89 CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_jacobian", &test_jacobian,
90 PETSC_NULLPTR);
91
92 CHKERR DMRegister_MoFEM("DMMOFEM");
93
94 Simple *si = m_field.getInterface<MoFEM::Simple>();
95
96 CHKERR si->getOptions();
97 CHKERR si->loadFile();
99 const int order = 2;
100 CHKERR si->setFieldOrder("x", order);
101
102 if (ale == PETSC_TRUE) {
104 CHKERR si->setFieldOrder("X", 2);
105 }
106 CHKERR si->setUp();
107
108 // create DM
109 DM dm;
110 CHKERR si->getDM(&dm);
111
112 // Projection on "x" field
113 {
114 Projection10NodeCoordsOnField ent_method(m_field, "x");
115 CHKERR m_field.loop_dofs("x", ent_method);
116 // CHKERR m_field.getInterface<FieldBlas>()->fieldScale(1.5, "x");
117 }
118
119 // Project coordinates on "X" field
120 if (ale == PETSC_TRUE) {
121 Projection10NodeCoordsOnField ent_method(m_field, "X");
122 CHKERR m_field.loop_dofs("X", ent_method);
123 // CHKERR m_field.getInterface<FieldBlas>()->fieldScale(1.5, "X");
124 }
125
126 boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
128 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
130 struct VolRule {
131 int operator()(int, int, int) const { return 2 * (order - 1); }
132 };
133 fe_lhs_ptr->getRuleHook = VolRule();
134 fe_rhs_ptr->getRuleHook = VolRule();
135
136 CHKERR DMMoFEMSNESSetJacobian(dm, si->getDomainFEName(), fe_lhs_ptr,
137 nullptr, nullptr);
138 CHKERR DMMoFEMSNESSetFunction(dm, si->getDomainFEName(), fe_rhs_ptr,
139 nullptr, nullptr);
140
142 boost::shared_ptr<map<int, BlockData>> block_sets_ptr =
143 boost::make_shared<map<int, BlockData>>();
144 (*block_sets_ptr)[0].iD = 0;
145 (*block_sets_ptr)[0].E = 1;
146 (*block_sets_ptr)[0].PoissonRatio = 0.25;
148 si->getDomainFEName(), 3, (*block_sets_ptr)[0].tEts);
149
150 // Parameters for density
151 const double rho_n = 2.0;
152 const double rho_0 = 0.5;
153
154 auto my_operators = [&](boost::shared_ptr<ForcesAndSourcesCore> &fe_lhs_ptr,
155 boost::shared_ptr<ForcesAndSourcesCore> &fe_rhs_ptr,
156 boost::shared_ptr<map<int, BlockData>>
157 &block_sets_ptr,
158 const std::string x_field,
159 const std::string X_field, const bool ale,
160 const bool field_disp) {
162
163 boost::shared_ptr<HookeElement::DataAtIntegrationPts> data_at_pts(
164 new HookeElement::DataAtIntegrationPts());
165 boost::shared_ptr<MatrixDouble> mat_coords_ptr =
166 boost::make_shared<MatrixDouble>();
167 boost::shared_ptr<VectorDouble> rho_at_gauss_pts_ptr =
168 boost::make_shared<VectorDouble>();
169 boost::shared_ptr<MatrixDouble> rho_grad_at_gauss_pts_ptr =
170 boost::make_shared<MatrixDouble>();
171
172 if (fe_lhs_ptr) {
173 if (ale == PETSC_FALSE) {
174 fe_lhs_ptr->getOpPtrVector().push_back(
175 new OpCalculateVectorFieldValues<3>(x_field, mat_coords_ptr));
176 fe_lhs_ptr->getOpPtrVector().push_back(new OpGetDensityField<false>(
177 x_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
178 rho_grad_at_gauss_pts_ptr));
179 fe_lhs_ptr->getOpPtrVector().push_back(
180 new HookeElement::OpCalculateStiffnessScaledByDensityField(
181 x_field, x_field, block_sets_ptr, data_at_pts,
182 rho_at_gauss_pts_ptr, rho_n, rho_0));
183 fe_lhs_ptr->getOpPtrVector().push_back(
184 new HookeElement::OpLhs_dx_dx<1>(x_field, x_field, data_at_pts));
185 } else {
186 fe_lhs_ptr->getOpPtrVector().push_back(
188 data_at_pts->HMat));
189 fe_lhs_ptr->getOpPtrVector().push_back(
190 new OpCalculateVectorFieldValues<3>(X_field, mat_coords_ptr));
191 fe_lhs_ptr->getOpPtrVector().push_back(new OpGetDensityField<true>(
192 X_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
193 rho_grad_at_gauss_pts_ptr));
194 fe_lhs_ptr->getOpPtrVector().push_back(
195 new HookeElement::OpCalculateStiffnessScaledByDensityField(
196 x_field, x_field, block_sets_ptr, data_at_pts,
197 rho_at_gauss_pts_ptr, rho_n, rho_0));
198 fe_lhs_ptr->getOpPtrVector().push_back(
200 data_at_pts->hMat));
201 fe_lhs_ptr->getOpPtrVector().push_back(
202 new HookeElement::OpCalculateStrainAle(x_field, x_field,
203 data_at_pts));
204 fe_lhs_ptr->getOpPtrVector().push_back(
205 new HookeElement::OpCalculateStress<1>(x_field, x_field,
206 data_at_pts)); // FIXME:
207 fe_lhs_ptr->getOpPtrVector().push_back(
208 new HookeElement::OpAleLhs_dx_dx<1>(x_field, x_field,
209 data_at_pts));
210 fe_lhs_ptr->getOpPtrVector().push_back(
211 new HookeElement::OpAleLhs_dx_dX<1>(x_field, X_field,
212 data_at_pts));
213 fe_lhs_ptr->getOpPtrVector().push_back(
214 new HookeElement::OpCalculateEnergy(X_field, X_field,
215 data_at_pts));
216 fe_lhs_ptr->getOpPtrVector().push_back(
217 new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
218 data_at_pts));
219 fe_lhs_ptr->getOpPtrVector().push_back(
220 new HookeElement::OpAleLhs_dX_dX<1>(X_field, X_field,
221 data_at_pts));
222 fe_lhs_ptr->getOpPtrVector().push_back(
223 new HookeElement::OpAleLhsPre_dX_dx<1>(X_field, x_field,
224 data_at_pts));
225 fe_lhs_ptr->getOpPtrVector().push_back(
226 new HookeElement::OpAleLhs_dX_dx(X_field, x_field, data_at_pts));
227 fe_lhs_ptr->getOpPtrVector().push_back(
228 new HookeElement::OpAleLhsWithDensity_dx_dX(
229 x_field, X_field, data_at_pts, rho_at_gauss_pts_ptr,
230 rho_grad_at_gauss_pts_ptr, rho_n, rho_0));
231 fe_lhs_ptr->getOpPtrVector().push_back(
232 new HookeElement::OpAleLhsWithDensity_dX_dX(
233 X_field, X_field, data_at_pts, rho_at_gauss_pts_ptr,
234 rho_grad_at_gauss_pts_ptr, rho_n, rho_0));
235 }
236 }
237
238 if (fe_rhs_ptr) {
239
240 if (ale == PETSC_FALSE) {
241 fe_rhs_ptr->getOpPtrVector().push_back(
243 data_at_pts->hMat));
244 fe_rhs_ptr->getOpPtrVector().push_back(
245 new OpCalculateVectorFieldValues<3>(x_field, mat_coords_ptr));
246 fe_rhs_ptr->getOpPtrVector().push_back(new OpGetDensityField<false>(
247 x_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
248 rho_grad_at_gauss_pts_ptr));
249 fe_rhs_ptr->getOpPtrVector().push_back(
250 new HookeElement::OpCalculateStiffnessScaledByDensityField(
251 x_field, x_field, block_sets_ptr, data_at_pts,
252 rho_at_gauss_pts_ptr, rho_n, rho_0));
253 if (field_disp) {
254 fe_rhs_ptr->getOpPtrVector().push_back(
255 new HookeElement::OpCalculateStrain<1>(x_field, x_field,
256 data_at_pts));
257 } else {
258 fe_rhs_ptr->getOpPtrVector().push_back(
259 new HookeElement::OpCalculateStrain<0>(x_field, x_field,
260 data_at_pts));
261 }
262 fe_rhs_ptr->getOpPtrVector().push_back(
263 new HookeElement::OpCalculateStress<1>(x_field, x_field,
264 data_at_pts));
265 fe_rhs_ptr->getOpPtrVector().push_back(
266 new HookeElement::OpRhs_dx(x_field, x_field, data_at_pts));
267 } else {
268 fe_rhs_ptr->getOpPtrVector().push_back(
270 data_at_pts->HMat));
271 fe_rhs_ptr->getOpPtrVector().push_back(
272 new OpCalculateVectorFieldValues<3>(X_field, mat_coords_ptr));
273 fe_rhs_ptr->getOpPtrVector().push_back(new OpGetDensityField<true>(
274 X_field, mat_coords_ptr, rho_at_gauss_pts_ptr,
275 rho_grad_at_gauss_pts_ptr));
276 fe_rhs_ptr->getOpPtrVector().push_back(
277 new HookeElement::OpCalculateStiffnessScaledByDensityField(
278 x_field, x_field, block_sets_ptr, data_at_pts,
279 rho_at_gauss_pts_ptr, rho_n, rho_0));
280 fe_rhs_ptr->getOpPtrVector().push_back(
282 data_at_pts->hMat));
283 fe_rhs_ptr->getOpPtrVector().push_back(
284 new HookeElement::OpCalculateStrainAle(x_field, x_field,
285 data_at_pts));
286 fe_rhs_ptr->getOpPtrVector().push_back(
287 new HookeElement::OpCalculateStress<1>(x_field, x_field,
288 data_at_pts));
289 fe_rhs_ptr->getOpPtrVector().push_back(
290 new HookeElement::OpAleRhs_dx(x_field, x_field, data_at_pts));
291 fe_rhs_ptr->getOpPtrVector().push_back(
292 new HookeElement::OpCalculateEnergy(X_field, X_field,
293 data_at_pts));
294 fe_rhs_ptr->getOpPtrVector().push_back(
295 new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
296 data_at_pts));
297 fe_rhs_ptr->getOpPtrVector().push_back(
298 new HookeElement::OpAleRhs_dX(X_field, X_field, data_at_pts));
299 }
300 }
302 };
303 CHKERR my_operators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr, "x", "X", ale,
304 false);
305 Vec x, f;
306 CHKERR DMCreateGlobalVector(dm, &x);
307 CHKERR VecDuplicate(x, &f);
308 CHKERR DMoFEMMeshToLocalVector(dm, x, INSERT_VALUES, SCATTER_FORWARD);
309
310 // CHKERR VecDuplicate(x, &dx);
311 // PetscRandom rctx;
312 // PetscRandomCreate(PETSC_COMM_WORLD, &rctx);
313 // VecSetRandom(dx, rctx);
314 // PetscRandomDestroy(&rctx);
315 // CHKERR DMoFEMMeshToGlobalVector(dm, x, INSERT_VALUES, SCATTER_REVERSE);
316
317 Mat A, fdA;
318 CHKERR DMCreateMatrix(dm, &A);
319 CHKERR MatDuplicate(A, MAT_DO_NOT_COPY_VALUES, &fdA);
320
321 if (test_jacobian == PETSC_TRUE) {
322 char testing_options[] =
323 "-snes_test_jacobian -snes_test_jacobian_display "
324 "-snes_no_convergence_test -snes_atol 0 -snes_rtol 0 -snes_max_it 1 "
325 "-pc_type none";
326 CHKERR PetscOptionsInsertString(NULL, testing_options);
327 } else {
328 char testing_options[] = "-snes_no_convergence_test -snes_atol 0 "
329 "-snes_rtol 0 -snes_max_it 1 -pc_type none";
330 CHKERR PetscOptionsInsertString(NULL, testing_options);
331 }
332
333 SNES snes;
334 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
335 MoFEM::SnesCtx *snes_ctx;
336 CHKERR DMMoFEMGetSnesCtx(dm, &snes_ctx);
337 CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
338 CHKERR SNESSetJacobian(snes, A, A, SnesMat, snes_ctx);
339 CHKERR SNESSetFromOptions(snes);
340
341 CHKERR SNESSolve(snes, NULL, x);
342
343 if (test_jacobian == PETSC_FALSE) {
344 double nrm_A0;
345 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A0);
346
347 char testing_options_fd[] = "-snes_fd";
348 CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
349
350 CHKERR SNESSetFunction(snes, f, SnesRhs, snes_ctx);
351 CHKERR SNESSetJacobian(snes, fdA, fdA, SnesMat, snes_ctx);
352 CHKERR SNESSetFromOptions(snes);
353
354 CHKERR SNESSolve(snes, NULL, x);
355 CHKERR MatAXPY(A, -1, fdA, SUBSET_NONZERO_PATTERN);
356
357 double nrm_A;
358 CHKERR MatNorm(A, NORM_INFINITY, &nrm_A);
359 PetscPrintf(PETSC_COMM_WORLD, "Matrix norms %3.4e %3.4e\n", nrm_A,
360 nrm_A / nrm_A0);
361 nrm_A /= nrm_A0;
362
363 const double tol = 1e-5;
364 if (nrm_A > tol) {
365 SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
366 "Difference between hand-calculated tangent matrix and finite "
367 "difference matrix is too big");
368 }
369 }
370
371 CHKERR VecDestroy(&x);
372 CHKERR VecDestroy(&f);
373 CHKERR MatDestroy(&A);
374 CHKERR MatDestroy(&fdA);
375 CHKERR SNESDestroy(&snes);
376
377 // destroy DM
378 CHKERR DMDestroy(&dm);
379 }
381
382 // finish work cleaning memory, getting statistics, etc
384
385 return 0;
386}
int main()
#define CATCH_ERRORS
Catch errors.
@ 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()
@ 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 ...
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
#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 ...
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:708
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:749
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition DMMoFEM.cpp:1084
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.
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:491
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:227
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
auto getFTensor1FromMat(M &data, int rr=0, int cc=0)
Get tensor rank 1 (vector) form data matrix.
static auto getFTensor0FromVec(V &data)
Get tensor rank 0 (scalar) form data vector.
constexpr AssemblyType A
Core (interface) class.
Definition Core.hpp:83
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition Core.cpp:68
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
Deprecated interface functions.
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Specialization for MatrixDouble vector field values calculation.
Projection of edge entities with one mid-node on hierarchical basis.
Simple interface for fast problem set-up.
Definition Simple.hpp:27
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:261
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:799
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:735
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:429
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:15
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference 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)
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 to volume elements.
int operator()(int, int, int) const