v0.16.0
Loading...
Searching...
No Matches
Classes | Functions | Variables
testing_jacobian_of_hook_scaled_with_density_element.cpp File Reference
#include <BasicFiniteElements.hpp>

Go to the source code of this file.

Classes

struct  OpGetDensityField< ALE >
 

Functions

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

Variables

static char help [] = "\n"
 

Function Documentation

◆ main()

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

Definition at line 72 of file testing_jacobian_of_hook_scaled_with_density_element.cpp.

72 {
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();
98 CHKERR si->addDomainField("x", H1, AINSWORTH_LEGENDRE_BASE, 3);
99 const int order = 2;
100 CHKERR si->setFieldOrder("x", order);
101
102 if (ale == PETSC_TRUE) {
103 CHKERR si->addDomainField("X", H1, AINSWORTH_LEGENDRE_BASE, 3);
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(
127 new VolumeElementForcesAndSourcesCore(m_field));
128 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
129 new VolumeElementForcesAndSourcesCore(m_field));
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
141 using BlockData = NonlinearElasticElement::BlockData;
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;
147 CHKERR m_field.get_finite_element_entities_by_dimension(
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) {
161 MoFEMFunctionBeginHot;
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(
187 new OpCalculateVectorFieldGradient<3, 3>(X_field,
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(
199 new OpCalculateVectorFieldGradient<3, 3>(x_field,
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(
242 new OpCalculateVectorFieldGradient<3, 3>(x_field,
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(
269 new OpCalculateVectorFieldGradient<3, 3>(X_field,
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(
281 new OpCalculateVectorFieldGradient<3, 3>(x_field,
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 }
301 MoFEMFunctionReturnHot(0);
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 }
380 CATCH_ERRORS;
381
382 // finish work cleaning memory, getting statistics, etc
383 MoFEM::Core::Finalize();
384
385 return 0;
386}
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition definitions.h:385
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
H1
@ H1
continuous field
Definition definitions.h:85
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
order
constexpr int order
Definition dg_projection.cpp:38
MoFEM::DMMoFEMSNESSetFunction
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
MoFEM::DMMoFEMSNESSetJacobian
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
MoFEM::DMoFEMMeshToLocalVector
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
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition DMMoFEM.cpp:43
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition DMMoFEM.cpp:1084
MoFEM::CoreInterface::get_finite_element_entities_by_dimension
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
MoFEM::CoreInterface::loop_dofs
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.
tol
double tol
Definition mesh_smoothing.cpp:23
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition MatrixFunction.hpp:64
MoFEM::SnesMat
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
MoFEM::SnesRhs
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
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition DeprecatedPetsc.hpp:182
A
constexpr AssemblyType A
Definition operators_tests.cpp:30
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:83
MoFEM::CoreTmp< 0 >::Initialize
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
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition Core.cpp:123
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Definition UserDataOperators.hpp:1711
MoFEM::OpCalculateVectorFieldValues
Specialization for MatrixDouble vector field values calculation.
Definition UserDataOperators.hpp:589
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition Projection10NodeCoordsOnField.hpp:24
MoFEM::Simple
Simple interface for fast problem set-up.
Definition Simple.hpp:27
MoFEM::Simple::addDomainField
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
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:799
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:735
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:429
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition SnesCtx.hpp:15
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
NonlinearElasticElement::BlockData
data for calculation heat conductivity and heat capacity elements
Definition HookeElement.hpp:32
VolRule
Set integration rule to volume elements.
Definition simple_interface.cpp:88
VolRule::operator()
int operator()(int, int, int) const
Definition simple_interface.cpp:89
help
static char help[]
Definition testing_jacobian_of_hook_scaled_with_density_element.cpp:16

Variable Documentation

◆ help

char help[] = "\n"
static

Definition at line 16 of file testing_jacobian_of_hook_scaled_with_density_element.cpp.

Generated by Doxygen 1.9.8 and hosted at University of Glasgow