1 #include <ceed.h> 2 #include <petsc.h> 3 #include "../problems/mooney-rivlin.h" 4 5 // Build libCEED context object 6 PetscErrorCode PhysicsContext_MR(MPI_Comm comm, Ceed ceed, Units *units, 7 CeedQFunctionContext *ctx) { 8 PetscErrorCode ierr; 9 Physics_MR phys; 10 11 PetscFunctionBegin; 12 13 ierr = PetscMalloc1(1, units); CHKERRQ(ierr); 14 ierr = PetscMalloc1(1, &phys); CHKERRQ(ierr); 15 ierr = ProcessPhysics_MR(comm, phys, *units); CHKERRQ(ierr); 16 CeedQFunctionContextCreate(ceed, ctx); 17 CeedQFunctionContextSetData(*ctx, CEED_MEM_HOST, CEED_COPY_VALUES, 18 sizeof(*phys), phys); 19 ierr = PetscFree(phys); CHKERRQ(ierr); 20 21 PetscFunctionReturn(0); 22 } 23 24 // Build libCEED smoother context object 25 PetscErrorCode PhysicsSmootherContext_MR(MPI_Comm comm, Ceed ceed, 26 CeedQFunctionContext ctx, CeedQFunctionContext *ctx_smoother) { 27 PetscErrorCode ierr; 28 PetscScalar nu_smoother = 0; 29 PetscBool nu_flag = PETSC_FALSE; 30 Physics_MR phys, phys_smoother; 31 32 PetscFunctionBegin; 33 34 ierr = PetscOptionsBegin(comm, NULL, 35 "Mooney Rivlin physical parameters for smoother", NULL); 36 CHKERRQ(ierr); 37 38 ierr = PetscOptionsScalar("-nu_smoother", "Poisson's ratio for smoother", 39 NULL, nu_smoother, &nu_smoother, &nu_flag); 40 CHKERRQ(ierr); 41 if (nu_smoother < 0 || 42 nu_smoother >= 0.5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, 43 "Mooney-Rivlin model requires Poisson ratio -nu option in [0, .5)"); 44 45 ierr = PetscOptionsEnd(); CHKERRQ(ierr); // End of setting Physics 46 47 if (nu_flag) { 48 // Copy context 49 CeedQFunctionContextGetData(ctx, CEED_MEM_HOST, &phys); 50 ierr = PetscMalloc1(1, &phys_smoother); CHKERRQ(ierr); 51 ierr = PetscMemcpy(phys_smoother, phys, sizeof(*phys)); CHKERRQ(ierr); 52 CeedQFunctionContextRestoreData(ctx, &phys); 53 // Create smoother context 54 CeedQFunctionContextCreate(ceed, ctx_smoother); 55 phys_smoother->lambda = 2 * (phys_smoother->mu_1 + phys_smoother->mu_2) * 56 nu_smoother / (1 - 2*nu_smoother); 57 CeedQFunctionContextSetData(*ctx_smoother, CEED_MEM_HOST, CEED_COPY_VALUES, 58 sizeof(*phys_smoother), phys_smoother); 59 ierr = PetscFree(phys_smoother); CHKERRQ(ierr); 60 } else { 61 *ctx_smoother = NULL; 62 } 63 64 PetscFunctionReturn(0); 65 } 66 67 // Process physics options - Mooney-Rivlin 68 PetscErrorCode ProcessPhysics_MR(MPI_Comm comm, Physics_MR phys, Units units) { 69 PetscErrorCode ierr; 70 PetscReal nu = -1; 71 phys->mu_1 = -1; 72 phys->mu_2 = -1; 73 phys->lambda = -1; 74 units->meter = 1; // 1 meter in scaled length units 75 units->second = 1; // 1 second in scaled time units 76 units->kilogram = 1; // 1 kilogram in scaled mass units 77 78 PetscFunctionBeginUser; 79 80 ierr = PetscOptionsBegin(comm, NULL, "Mooney Rivlin physical parameters", NULL); 81 CHKERRQ(ierr); 82 83 ierr = PetscOptionsScalar("-mu_1", "Material Property mu_1", NULL, 84 phys->mu_1, &phys->mu_1, NULL); CHKERRQ(ierr); 85 if (phys->mu_1 < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, 86 "Mooney-Rivlin model requires non-negative -mu_1 option (Pa)"); 87 88 ierr = PetscOptionsScalar("-mu_2", "Material Property mu_2", NULL, 89 phys->mu_2, &phys->mu_2, NULL); CHKERRQ(ierr); 90 if (phys->mu_2 < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, 91 "Mooney-Rivlin model requires non-negative -mu_2 option (Pa)"); 92 93 ierr = PetscOptionsScalar("-nu", "Poisson ratio", NULL, 94 nu, &nu, NULL); CHKERRQ(ierr); 95 if (nu < 0 || nu >= 0.5) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, 96 "Mooney-Rivlin model requires Poisson ratio -nu option in [0, .5)"); 97 phys->lambda = 2 * (phys->mu_1 + phys->mu_2) * nu / (1 - 2*nu); 98 99 ierr = PetscOptionsScalar("-units_meter", "1 meter in scaled length units", 100 NULL, units->meter, &units->meter, NULL); 101 CHKERRQ(ierr); 102 units->meter = fabs(units->meter); 103 104 ierr = PetscOptionsScalar("-units_second", "1 second in scaled time units", 105 NULL, units->second, &units->second, NULL); 106 CHKERRQ(ierr); 107 units->second = fabs(units->second); 108 109 ierr = PetscOptionsScalar("-units_kilogram", "1 kilogram in scaled mass units", 110 NULL, units->kilogram, &units->kilogram, NULL); 111 CHKERRQ(ierr); 112 units->kilogram = fabs(units->kilogram); 113 114 ierr = PetscOptionsEnd(); CHKERRQ(ierr); // End of setting Physics 115 116 // Define derived units 117 units->Pascal = units->kilogram / (units->meter * PetscSqr(units->second)); 118 119 // Scale material parameters based on units of Pa 120 phys->mu_1 *= units->Pascal; 121 phys->mu_2 *= units->Pascal; 122 phys->lambda *= units->Pascal; 123 124 PetscFunctionReturn(0); 125 };