v0.14.0
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 
4 Testing implementation of Hook element by verifying tangent stiffness matrix.
5 Test like this is an example of how to verify the implementation of Jacobian.
6 
7 */
8 
9 
10 
11 #include <BasicFiniteElements.hpp>
12 
13 using namespace boost::numeric;
14 using namespace MoFEM;
15 
16 static char help[] = "\n";
17 
18 template <bool ALE>
19 struct 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) {
34  MoFEMFunctionBegin;
35  if (row_type != MBVERTEX)
36  MoFEMFunctionReturnHot(0);
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 
44  auto t_rho = getFTensor0FromVec(*rhoAtGaussPtsPtr);
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 
53  FTensor::Index<'i', 3> i;
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 
66  MoFEMFunctionReturn(0);
67  }
68 };
69 
70 int 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();
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 }
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
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.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition: Core.hpp:82
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
MoFEM::OpCalculateVectorFieldValues
Get values at integration pts for tensor filed rank 1, i.e. vector field.
Definition: UserDataOperators.hpp:466
MoFEM::Types::MatrixDouble
UBlasMatrix< double > MatrixDouble
Definition: Types.hpp:77
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
A
constexpr AssemblyType A
Definition: operators_tests.cpp:30
VolRule
Set integration rule to volume elements.
Definition: simple_interface.cpp:88
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:523
BasicFiniteElements.hpp
MoFEM::Projection10NodeCoordsOnField
Projection of edge entities with one mid-node on hierarchical basis.
Definition: Projection10NodeCoordsOnField.hpp:24
MoFEM::CoreTmp< 0 >::Finalize
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition: DeprecatedCoreInterface.hpp:16
MoFEM::Interface
DeprecatedCoreInterface Interface
Definition: Interface.hpp:2010
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:747
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:10
OpGetDensityField::OpGetDensityField
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)
Definition: testing_jacobian_of_hook_scaled_with_density_element.cpp:24
MoFEM::getFTensor0FromVec
static auto getFTensor0FromVec(ublas::vector< T, A > &data)
Get tensor rank 0 (scalar) form data vector.
Definition: Templates.hpp:135
HookeElement
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::DMMoFEMGetSnesCtx
PetscErrorCode DMMoFEMGetSnesCtx(DM dm, MoFEM::SnesCtx **snes_ctx)
get MoFEM::SnesCtx data structure
Definition: DMMoFEM.cpp:1094
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:368
OpGetDensityField::doWork
MoFEMErrorCode doWork(int row_side, EntityType row_type, HookeElement::EntData &row_data)
Definition: testing_jacobian_of_hook_scaled_with_density_element.cpp:32
VolUserDataOperator
VolumeElementForcesAndSourcesCore::UserDataOperator VolUserDataOperator
Definition: elasticity.cpp:76
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:759
main
int main(int argc, char *argv[])
Definition: testing_jacobian_of_hook_scaled_with_density_element.cpp:70
OpGetDensityField::rhoGradAtGaussPtsPtr
boost::shared_ptr< MatrixDouble > rhoGradAtGaussPtsPtr
Definition: testing_jacobian_of_hook_scaled_with_density_element.cpp:23
help
static char help[]
Definition: testing_jacobian_of_hook_scaled_with_density_element.cpp:16
MoFEM::VolumeElementForcesAndSourcesCore
Volume finite element base.
Definition: VolumeElementForcesAndSourcesCore.hpp:26
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator
Definition: VolumeElementForcesAndSourcesCore.hpp:108
EntData
EntitiesFieldData::EntData EntData
Definition: child_and_parent.cpp:37
FTensor::Index< 'i', 3 >
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:545
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::OpCalculateVectorFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1535
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::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
NonlinearElasticElement::BlockData::iD
int iD
Definition: HookeElement.hpp:33
CATCH_ERRORS
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
MoFEM::Core
CoreTmp< 0 > Core
Definition: Core.hpp:1148
HenckyOps::f
auto f
Definition: HenckyOps.hpp:15
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
OpGetDensityField
Definition: testing_jacobian_of_hook_scaled_with_density_element.cpp:19
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition: MatrixFunction.hpp:66
OpGetDensityField::rhoAtGaussPtsPtr
boost::shared_ptr< VectorDouble > rhoAtGaussPtsPtr
Definition: testing_jacobian_of_hook_scaled_with_density_element.cpp:22
MoFEM::BlockData
Definition: MeshsetsManager.cpp:755
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
MoFEM::SnesCtx
Interface for nonlinear (SNES) solver.
Definition: SnesCtx.hpp:13
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:139
OpGetDensityField::matCoordsPtr
boost::shared_ptr< MatrixDouble > matCoordsPtr
Definition: testing_jacobian_of_hook_scaled_with_density_element.cpp:21
MoFEM::Simple::setUp
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:683
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
BlockData
NonlinearElasticElement::BlockData BlockData
Definition: elasticity.cpp:74
tol
double tol
Definition: mesh_smoothing.cpp:26
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:718
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
Generated by Doxygen 1.8.17 and hosted at University of Glasgow