72int main(
int argc,
char *argv[]) {
78 moab::Core mb_instance;
79 moab::Interface &moab = mb_instance;
86 PetscBool ale = PETSC_FALSE;
88 PetscBool test_jacobian = PETSC_FALSE;
102 if (ale == PETSC_TRUE) {
120 if (ale == PETSC_TRUE) {
126 boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
128 boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
133 fe_lhs_ptr->getRuleHook =
VolRule();
134 fe_rhs_ptr->getRuleHook =
VolRule();
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;
151 const double rho_n = 2.0;
152 const double rho_0 = 0.5;
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>>
158 const std::string x_field,
159 const std::string X_field,
const bool ale,
160 const bool field_disp) {
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>();
173 if (ale == PETSC_FALSE) {
174 fe_lhs_ptr->getOpPtrVector().push_back(
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));
186 fe_lhs_ptr->getOpPtrVector().push_back(
189 fe_lhs_ptr->getOpPtrVector().push_back(
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(
201 fe_lhs_ptr->getOpPtrVector().push_back(
202 new HookeElement::OpCalculateStrainAle(x_field, x_field,
204 fe_lhs_ptr->getOpPtrVector().push_back(
205 new HookeElement::OpCalculateStress<1>(x_field, x_field,
207 fe_lhs_ptr->getOpPtrVector().push_back(
208 new HookeElement::OpAleLhs_dx_dx<1>(x_field, x_field,
210 fe_lhs_ptr->getOpPtrVector().push_back(
211 new HookeElement::OpAleLhs_dx_dX<1>(x_field, X_field,
213 fe_lhs_ptr->getOpPtrVector().push_back(
214 new HookeElement::OpCalculateEnergy(X_field, X_field,
216 fe_lhs_ptr->getOpPtrVector().push_back(
217 new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
219 fe_lhs_ptr->getOpPtrVector().push_back(
220 new HookeElement::OpAleLhs_dX_dX<1>(X_field, X_field,
222 fe_lhs_ptr->getOpPtrVector().push_back(
223 new HookeElement::OpAleLhsPre_dX_dx<1>(X_field, x_field,
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));
240 if (ale == PETSC_FALSE) {
241 fe_rhs_ptr->getOpPtrVector().push_back(
244 fe_rhs_ptr->getOpPtrVector().push_back(
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));
254 fe_rhs_ptr->getOpPtrVector().push_back(
255 new HookeElement::OpCalculateStrain<1>(x_field, x_field,
258 fe_rhs_ptr->getOpPtrVector().push_back(
259 new HookeElement::OpCalculateStrain<0>(x_field, x_field,
262 fe_rhs_ptr->getOpPtrVector().push_back(
263 new HookeElement::OpCalculateStress<1>(x_field, x_field,
265 fe_rhs_ptr->getOpPtrVector().push_back(
266 new HookeElement::OpRhs_dx(x_field, x_field, data_at_pts));
268 fe_rhs_ptr->getOpPtrVector().push_back(
271 fe_rhs_ptr->getOpPtrVector().push_back(
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(
283 fe_rhs_ptr->getOpPtrVector().push_back(
284 new HookeElement::OpCalculateStrainAle(x_field, x_field,
286 fe_rhs_ptr->getOpPtrVector().push_back(
287 new HookeElement::OpCalculateStress<1>(x_field, x_field,
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,
294 fe_rhs_ptr->getOpPtrVector().push_back(
295 new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
297 fe_rhs_ptr->getOpPtrVector().push_back(
298 new HookeElement::OpAleRhs_dX(X_field, X_field, data_at_pts));
303 CHKERR my_operators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr,
"x",
"X", ale,
306 CHKERR DMCreateGlobalVector(dm, &x);
307 CHKERR VecDuplicate(x, &f);
319 CHKERR MatDuplicate(
A, MAT_DO_NOT_COPY_VALUES, &fdA);
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 "
326 CHKERR PetscOptionsInsertString(NULL, testing_options);
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);
334 CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
339 CHKERR SNESSetFromOptions(snes);
341 CHKERR SNESSolve(snes, NULL, x);
343 if (test_jacobian == PETSC_FALSE) {
345 CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A0);
347 char testing_options_fd[] =
"-snes_fd";
348 CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
352 CHKERR SNESSetFromOptions(snes);
354 CHKERR SNESSolve(snes, NULL, x);
355 CHKERR MatAXPY(
A, -1, fdA, SUBSET_NONZERO_PATTERN);
358 CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A);
359 PetscPrintf(PETSC_COMM_WORLD,
"Matrix norms %3.4e %3.4e\n", nrm_A,
363 const double tol = 1e-5;
366 "Difference between hand-calculated tangent matrix and finite "
367 "difference matrix is too big");
375 CHKERR SNESDestroy(&snes);