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