1 // Copyright (c) 2017-2026, 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
PhysicsContext_MR(MPI_Comm comm,Ceed ceed,Units * units,CeedQFunctionContext * ctx)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
PhysicsSmootherContext_MR(MPI_Comm comm,Ceed ceed,CeedQFunctionContext ctx,CeedQFunctionContext * ctx_smoother)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
ProcessPhysics_MR(MPI_Comm comm,Physics_MR phys,Units units)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 };