v0.14.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 70 of file testing_jacobian_of_hook_scaled_with_density_element.cpp.

70 {
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();
96 CHKERR si->addDomainField("x", H1, AINSWORTH_LEGENDRE_BASE, 3);
97 const int order = 2;
98 CHKERR si->setFieldOrder("x", order);
99
100 if (ale == PETSC_TRUE) {
101 CHKERR si->addDomainField("X", H1, AINSWORTH_LEGENDRE_BASE, 3);
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(
125 new VolumeElementForcesAndSourcesCore(m_field));
126 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
127 new VolumeElementForcesAndSourcesCore(m_field));
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
139 using BlockData = NonlinearElasticElement::BlockData;
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;
145 CHKERR m_field.get_finite_element_entities_by_dimension(
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) {
159 MoFEMFunctionBeginHot;
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(
185 new OpCalculateVectorFieldGradient<3, 3>(X_field,
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(
197 new OpCalculateVectorFieldGradient<3, 3>(x_field,
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(
240 new OpCalculateVectorFieldGradient<3, 3>(x_field,
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(
267 new OpCalculateVectorFieldGradient<3, 3>(X_field,
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(
279 new OpCalculateVectorFieldGradient<3, 3>(x_field,
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 }
299 MoFEMFunctionReturnHot(0);
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 }
378 CATCH_ERRORS;
379
380 // finish work cleaning memory, getting statistics, etc
381 MoFEM::Core::Finalize();
382
383 return 0;
384}
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
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:447
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:535
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
order
constexpr int order
Definition: dg_projection.cpp:18
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:704
MoFEM::DMoFEMMeshToLocalVector
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
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:745
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:47
MoFEM::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1080
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:27
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
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:136
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:27
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:82
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:72
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1520
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:450
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:264
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:667
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:473
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:611
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:327
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce 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.
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.5 and hosted at University of Glasgow