70int main(
int argc, 
char *argv[]) {
 
   76  moab::Core mb_instance;              
 
   77  moab::Interface &moab = mb_instance; 
 
   84    PetscBool ale = PETSC_FALSE;
 
   86    PetscBool test_jacobian = PETSC_FALSE;
 
  100    if (ale == PETSC_TRUE) {
 
  118    if (ale == PETSC_TRUE) {
 
  124    boost::shared_ptr<ForcesAndSourcesCore> fe_lhs_ptr(
 
  126    boost::shared_ptr<ForcesAndSourcesCore> fe_rhs_ptr(
 
  131    fe_lhs_ptr->getRuleHook = 
VolRule();
 
  132    fe_rhs_ptr->getRuleHook = 
VolRule();
 
  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;
 
  149    const double rho_n = 2.0;
 
  150    const double rho_0 = 0.5;
 
  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>>
 
  156                            const std::string x_field,
 
  157                            const std::string X_field, 
const bool ale,
 
  158                            const bool field_disp) {
 
  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>();
 
  171        if (ale == PETSC_FALSE) {
 
  172          fe_lhs_ptr->getOpPtrVector().push_back(
 
  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));
 
  184          fe_lhs_ptr->getOpPtrVector().push_back(
 
  187          fe_lhs_ptr->getOpPtrVector().push_back(
 
  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(
 
  199          fe_lhs_ptr->getOpPtrVector().push_back(
 
  200              new HookeElement::OpCalculateStrainAle(x_field, x_field,
 
  202          fe_lhs_ptr->getOpPtrVector().push_back(
 
  203              new HookeElement::OpCalculateStress<1>(x_field, x_field,
 
  205          fe_lhs_ptr->getOpPtrVector().push_back(
 
  206              new HookeElement::OpAleLhs_dx_dx<1>(x_field, x_field,
 
  208          fe_lhs_ptr->getOpPtrVector().push_back(
 
  209              new HookeElement::OpAleLhs_dx_dX<1>(x_field, X_field,
 
  211          fe_lhs_ptr->getOpPtrVector().push_back(
 
  212              new HookeElement::OpCalculateEnergy(X_field, X_field,
 
  214          fe_lhs_ptr->getOpPtrVector().push_back(
 
  215              new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
 
  217          fe_lhs_ptr->getOpPtrVector().push_back(
 
  218              new HookeElement::OpAleLhs_dX_dX<1>(X_field, X_field,
 
  220          fe_lhs_ptr->getOpPtrVector().push_back(
 
  221              new HookeElement::OpAleLhsPre_dX_dx<1>(X_field, x_field,
 
  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));
 
  238        if (ale == PETSC_FALSE) {
 
  239          fe_rhs_ptr->getOpPtrVector().push_back(
 
  242          fe_rhs_ptr->getOpPtrVector().push_back(
 
  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));
 
  252            fe_rhs_ptr->getOpPtrVector().push_back(
 
  253                new HookeElement::OpCalculateStrain<1>(x_field, x_field,
 
  256            fe_rhs_ptr->getOpPtrVector().push_back(
 
  257                new HookeElement::OpCalculateStrain<0>(x_field, x_field,
 
  260          fe_rhs_ptr->getOpPtrVector().push_back(
 
  261              new HookeElement::OpCalculateStress<1>(x_field, x_field,
 
  263          fe_rhs_ptr->getOpPtrVector().push_back(
 
  264              new HookeElement::OpRhs_dx(x_field, x_field, data_at_pts));
 
  266          fe_rhs_ptr->getOpPtrVector().push_back(
 
  269          fe_rhs_ptr->getOpPtrVector().push_back(
 
  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(
 
  281          fe_rhs_ptr->getOpPtrVector().push_back(
 
  282              new HookeElement::OpCalculateStrainAle(x_field, x_field,
 
  284          fe_rhs_ptr->getOpPtrVector().push_back(
 
  285              new HookeElement::OpCalculateStress<1>(x_field, x_field,
 
  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,
 
  292          fe_rhs_ptr->getOpPtrVector().push_back(
 
  293              new HookeElement::OpCalculateEshelbyStress(X_field, X_field,
 
  295          fe_rhs_ptr->getOpPtrVector().push_back(
 
  296              new HookeElement::OpAleRhs_dX(X_field, X_field, data_at_pts));
 
  301    CHKERR my_operators(fe_lhs_ptr, fe_rhs_ptr, block_sets_ptr, 
"x", 
"X", ale,
 
  304    CHKERR DMCreateGlobalVector(dm, &x);
 
  305    CHKERR VecDuplicate(x, &f);
 
  317    CHKERR MatDuplicate(
A, MAT_DO_NOT_COPY_VALUES, &fdA);
 
  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 " 
  324      CHKERR PetscOptionsInsertString(NULL, testing_options);
 
  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);
 
  332    CHKERR SNESCreate(PETSC_COMM_WORLD, &snes);
 
  337    CHKERR SNESSetFromOptions(snes);
 
  339    CHKERR SNESSolve(snes, NULL, x);
 
  341    if (test_jacobian == PETSC_FALSE) {
 
  343      CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A0);
 
  345      char testing_options_fd[] = 
"-snes_fd";
 
  346      CHKERR PetscOptionsInsertString(NULL, testing_options_fd);
 
  350      CHKERR SNESSetFromOptions(snes);
 
  352      CHKERR SNESSolve(snes, NULL, x);
 
  353      CHKERR MatAXPY(
A, -1, fdA, SUBSET_NONZERO_PATTERN);
 
  356      CHKERR MatNorm(
A, NORM_INFINITY, &nrm_A);
 
  357      PetscPrintf(PETSC_COMM_WORLD, 
"Matrix norms %3.4e %3.4e\n", nrm_A,
 
  361      const double tol = 1e-5;
 
  364                "Difference between hand-calculated tangent matrix and finite " 
  365                "difference matrix is too big");
 
  373    CHKERR SNESDestroy(&snes);