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