1 /// @file 2 /// Command line option processing for solid mechanics example using PETSc 3 4 #include "../include/cl-options.h" 5 6 // ----------------------------------------------------------------------------- 7 // Process command line options 8 // ----------------------------------------------------------------------------- 9 // Process general command line options 10 PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx) { 11 PetscErrorCode ierr; 12 PetscBool ceed_flag = PETSC_FALSE; 13 14 PetscFunctionBeginUser; 15 16 ierr = PetscOptionsBegin(comm, NULL, 17 "Elasticity / Hyperelasticity in PETSc with libCEED", 18 NULL); CHKERRQ(ierr); 19 20 ierr = PetscOptionsString("-ceed", "CEED resource specifier", 21 NULL, app_ctx->ceed_resource, app_ctx->ceed_resource, 22 sizeof(app_ctx->ceed_resource), &ceed_flag); 23 CHKERRQ(ierr); 24 25 ierr = PetscStrncpy(app_ctx->output_dir, ".", 2); 26 CHKERRQ(ierr); // Default - current directory 27 ierr = PetscOptionsString("-output_dir", "Output directory", 28 NULL, app_ctx->output_dir, app_ctx->output_dir, 29 sizeof(app_ctx->output_dir), NULL); CHKERRQ(ierr); 30 31 app_ctx->degree = 3; 32 ierr = PetscOptionsInt("-degree", "Polynomial degree of tensor product basis", 33 NULL, app_ctx->degree, &app_ctx->degree, NULL); 34 CHKERRQ(ierr); 35 36 app_ctx->q_extra = 0; 37 ierr = PetscOptionsInt("-qextra", "Number of extra quadrature points", 38 NULL, app_ctx->q_extra, &app_ctx->q_extra, NULL); 39 CHKERRQ(ierr); 40 41 ierr = PetscOptionsString("-mesh", "Read mesh from file", NULL, 42 app_ctx->mesh_file, app_ctx->mesh_file, 43 sizeof(app_ctx->mesh_file), NULL); CHKERRQ(ierr); 44 45 app_ctx->problem_choice = ELAS_LINEAR; // Default - Linear Elasticity 46 ierr = PetscOptionsEnum("-problem", 47 "Solves Elasticity & Hyperelasticity Problems", 48 NULL, problemTypes, (PetscEnum)app_ctx->problem_choice, 49 (PetscEnum *)&app_ctx->problem_choice, NULL); 50 CHKERRQ(ierr); 51 app_ctx->name = problemTypes[app_ctx->problem_choice]; 52 app_ctx->name_for_disp = problemTypesForDisp[app_ctx->problem_choice]; 53 54 app_ctx->num_increments = app_ctx->problem_choice == ELAS_LINEAR ? 1 : 10; 55 ierr = PetscOptionsInt("-num_steps", "Number of pseudo-time steps", 56 NULL, app_ctx->num_increments, &app_ctx->num_increments, 57 NULL); CHKERRQ(ierr); 58 59 app_ctx->forcing_choice = FORCE_NONE; // Default - no forcing term 60 ierr = PetscOptionsEnum("-forcing", "Set forcing function option", NULL, 61 forcing_types, (PetscEnum)app_ctx->forcing_choice, 62 (PetscEnum *)&app_ctx->forcing_choice, NULL); 63 CHKERRQ(ierr); 64 65 PetscInt max_n = 3; 66 app_ctx->forcing_vector[0] = 0; 67 app_ctx->forcing_vector[1] = -1; 68 app_ctx->forcing_vector[2] = 0; 69 ierr = PetscOptionsScalarArray("-forcing_vec", 70 "Direction to apply constant force", NULL, 71 app_ctx->forcing_vector, &max_n, NULL); 72 CHKERRQ(ierr); 73 74 if ((app_ctx->problem_choice == ELAS_FSInitial_NH1 || 75 app_ctx->problem_choice == ELAS_FSInitial_NH2 || 76 app_ctx->problem_choice == ELAS_FSCurrent_NH1 || 77 app_ctx->problem_choice == ELAS_FSCurrent_NH2 || 78 app_ctx->problem_choice == ELAS_FSInitial_MR1) && 79 app_ctx->forcing_choice == FORCE_CONST) 80 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, 81 "Cannot use constant forcing and finite strain formulation. " 82 "Constant forcing in reference frame currently unavaliable."); 83 84 // Dirichlet boundary conditions 85 app_ctx->bc_clamp_count = 16; 86 ierr = PetscOptionsIntArray("-bc_clamp", 87 "Face IDs to apply incremental Dirichlet BC", 88 NULL, app_ctx->bc_clamp_faces, &app_ctx->bc_clamp_count, 89 NULL); CHKERRQ(ierr); 90 // Set vector for each clamped BC 91 for (PetscInt i = 0; i < app_ctx->bc_clamp_count; i++) { 92 // Translation vector 93 char option_name[25]; 94 const size_t nclamp_params = sizeof(app_ctx->bc_clamp_max[0])/sizeof( 95 app_ctx->bc_clamp_max[0][0]); 96 for (PetscInt j = 0; j < nclamp_params; j++) 97 app_ctx->bc_clamp_max[i][j] = 0.; 98 99 snprintf(option_name, sizeof option_name, "-bc_clamp_%d_translate", 100 app_ctx->bc_clamp_faces[i]); 101 max_n = 3; 102 ierr = PetscOptionsScalarArray(option_name, 103 "Vector to translate clamped end by", NULL, 104 app_ctx->bc_clamp_max[i], &max_n, NULL); 105 CHKERRQ(ierr); 106 107 // Rotation vector 108 max_n = 5; 109 snprintf(option_name, sizeof option_name, "-bc_clamp_%d_rotate", 110 app_ctx->bc_clamp_faces[i]); 111 ierr = PetscOptionsScalarArray(option_name, 112 "Vector with axis of rotation and rotation, in radians", 113 NULL, &app_ctx->bc_clamp_max[i][3], &max_n, NULL); 114 CHKERRQ(ierr); 115 116 // Normalize 117 PetscScalar norm = sqrt(app_ctx->bc_clamp_max[i][3]*app_ctx->bc_clamp_max[i][3] 118 + app_ctx->bc_clamp_max[i][4]*app_ctx->bc_clamp_max[i][4] 119 + app_ctx->bc_clamp_max[i][5]*app_ctx->bc_clamp_max[i][5]); 120 if (fabs(norm) < 1e-16) 121 norm = 1; 122 for (PetscInt j = 0; j < 3; j++) 123 app_ctx->bc_clamp_max[i][3 + j] /= norm; 124 } 125 126 // Neumann boundary conditions 127 app_ctx->bc_traction_count = 16; 128 ierr = PetscOptionsIntArray("-bc_traction", 129 "Face IDs to apply traction (Neumann) BC", 130 NULL, app_ctx->bc_traction_faces, 131 &app_ctx->bc_traction_count, NULL); CHKERRQ(ierr); 132 // Set vector for each traction BC 133 for (PetscInt i = 0; i < app_ctx->bc_traction_count; i++) { 134 // Translation vector 135 char option_name[25]; 136 for (PetscInt j = 0; j < 3; j++) 137 app_ctx->bc_traction_vector[i][j] = 0.; 138 139 snprintf(option_name, sizeof option_name, "-bc_traction_%d", 140 app_ctx->bc_traction_faces[i]); 141 max_n = 3; 142 PetscBool set = false; 143 ierr = PetscOptionsScalarArray(option_name, 144 "Traction vector for constrained face", NULL, 145 app_ctx->bc_traction_vector[i], &max_n, &set); 146 CHKERRQ(ierr); 147 148 if (!set) 149 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, 150 "Traction vector must be set for all traction boundary conditions."); 151 } 152 153 app_ctx->multigrid_choice = MULTIGRID_LOGARITHMIC; 154 ierr = PetscOptionsEnum("-multigrid", "Set multigrid type option", NULL, 155 multigrid_types, (PetscEnum)app_ctx->multigrid_choice, 156 (PetscEnum *)&app_ctx->multigrid_choice, NULL); 157 CHKERRQ(ierr); 158 159 app_ctx->test_mode = PETSC_FALSE; 160 ierr = PetscOptionsBool("-test", 161 "Testing mode (do not print unless error is large)", 162 NULL, app_ctx->test_mode, &(app_ctx->test_mode), NULL); 163 CHKERRQ(ierr); 164 165 app_ctx->expect_final_strain = -1.; 166 ierr = PetscOptionsReal("-expect_final_strain_energy", 167 "Expect final strain energy close to this value.", 168 NULL, app_ctx->expect_final_strain, &app_ctx->expect_final_strain, NULL); 169 CHKERRQ(ierr); 170 171 app_ctx->test_tol = 1e-8; 172 ierr = PetscOptionsReal("-expect_final_state_rtol", 173 "Relative tolerance for final strain energy test", 174 NULL, app_ctx->test_tol, &app_ctx->test_tol, NULL); 175 CHKERRQ(ierr); 176 177 app_ctx->view_soln = PETSC_FALSE; 178 ierr = PetscOptionsBool("-view_soln", "Write out solution vector for viewing", 179 NULL, app_ctx->view_soln, &(app_ctx->view_soln), NULL); 180 CHKERRQ(ierr); 181 182 app_ctx->view_final_soln = PETSC_FALSE; 183 ierr = PetscOptionsBool("-view_final_soln", 184 "Write out final solution vector for viewing", 185 NULL, app_ctx->view_final_soln, &(app_ctx->view_final_soln), 186 NULL); CHKERRQ(ierr); 187 CHKERRQ(ierr); 188 189 PetscBool set; 190 char energy_viewer_filename[PETSC_MAX_PATH_LEN] = ""; 191 ierr = PetscOptionsString("-strain_energy_monitor", 192 "Print out current strain energy at every load increment", 193 NULL, energy_viewer_filename, 194 energy_viewer_filename, sizeof(energy_viewer_filename), 195 &set); CHKERRQ(ierr); 196 if (set) { 197 ierr = PetscViewerASCIIOpen(comm, energy_viewer_filename, 198 &app_ctx->energy_viewer); CHKERRQ(ierr); 199 ierr = PetscViewerASCIIPrintf(app_ctx->energy_viewer, "increment,energy\n"); 200 CHKERRQ(ierr); 201 // Initial configuration is base energy state; this may not be true if we extend in the future to 202 // initially loaded configurations (because a truly at-rest initial state may not be realizable). 203 ierr = PetscViewerASCIIPrintf(app_ctx->energy_viewer, "%f,%e\n", 0., 0.); 204 CHKERRQ(ierr); 205 } 206 ierr = PetscOptionsEnd(); CHKERRQ(ierr); // End of setting AppCtx 207 208 // Check for all required values set 209 if (app_ctx->test_mode) { 210 if (app_ctx->forcing_choice == FORCE_NONE && !app_ctx->bc_clamp_count) 211 app_ctx->forcing_choice = FORCE_MMS; 212 } 213 if (!app_ctx->bc_clamp_count && app_ctx->forcing_choice != FORCE_MMS) { 214 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "-boundary options needed"); 215 } 216 217 // Provide default ceed resource if not specified 218 if (!ceed_flag) { 219 const char *ceed_resource = "/cpu/self"; 220 strncpy(app_ctx->ceed_resource, ceed_resource, 10); 221 } 222 223 // Determine number of levels 224 switch (app_ctx->multigrid_choice) { 225 case MULTIGRID_LOGARITHMIC: 226 app_ctx->num_levels = ceil(log(app_ctx->degree)/log(2)) + 1; 227 break; 228 case MULTIGRID_UNIFORM: 229 app_ctx->num_levels = app_ctx->degree; 230 break; 231 case MULTIGRID_NONE: 232 app_ctx->num_levels = 1; 233 break; 234 } 235 236 // Populate array of degrees for each level for multigrid 237 ierr = PetscMalloc1(app_ctx->num_levels, &(app_ctx->level_degrees)); 238 CHKERRQ(ierr); 239 240 switch (app_ctx->multigrid_choice) { 241 case MULTIGRID_LOGARITHMIC: 242 for (int i = 0; i < app_ctx->num_levels-1; i++) 243 app_ctx->level_degrees[i] = pow(2,i); 244 app_ctx->level_degrees[app_ctx->num_levels-1] = app_ctx->degree; 245 break; 246 case MULTIGRID_UNIFORM: 247 for (int i = 0; i < app_ctx->num_levels; i++) 248 app_ctx->level_degrees[i] = i + 1; 249 break; 250 case MULTIGRID_NONE: 251 app_ctx->level_degrees[0] = app_ctx->degree; 252 break; 253 } 254 255 PetscFunctionReturn(0); 256 }; 257