1 /// @file 2 /// Boundary condition functions for solid mechanics example using PETSc 3 4 #include "../include/boundary.h" 5 6 // ----------------------------------------------------------------------------- 7 // Boundary Functions 8 // ----------------------------------------------------------------------------- 9 // Note: If additional boundary conditions are added, an update is needed in 10 // elasticity.h for the boundaryOptions variable. 11 12 // BCMMS - boundary function 13 // Values on all points of the mesh is set based on given solution below 14 // for u[0], u[1], u[2] 15 PetscErrorCode BCMMS(PetscInt dim, PetscReal load_increment, 16 const PetscReal coords[], PetscInt num_comp_u, 17 PetscScalar *u, void *ctx) { 18 PetscScalar x = coords[0]; 19 PetscScalar y = coords[1]; 20 PetscScalar z = coords[2]; 21 22 PetscFunctionBeginUser; 23 24 u[0] = exp(2*x)*sin(3*y)*cos(4*z) / 1e8 * load_increment; 25 u[1] = exp(3*y)*sin(4*z)*cos(2*x) / 1e8 * load_increment; 26 u[2] = exp(4*z)*sin(2*x)*cos(3*y) / 1e8 * load_increment; 27 28 PetscFunctionReturn(0); 29 }; 30 31 #ifndef M_PI 32 # define M_PI 3.14159265358979323846 33 #endif 34 35 // BCClamp - fix boundary values with affine transformation at fraction of load 36 // increment 37 PetscErrorCode BCClamp(PetscInt dim, PetscReal load_increment, 38 const PetscReal coords[], PetscInt num_comp_u, 39 PetscScalar *u, void *ctx) { 40 PetscScalar x = coords[0]; 41 PetscScalar y = coords[1]; 42 PetscScalar z = coords[2]; 43 PetscScalar (*clampMax) = (PetscScalar(*))ctx; 44 45 PetscFunctionBeginUser; 46 PetscScalar 47 // Translation vector 48 lx = clampMax[0]*load_increment, 49 ly = clampMax[1]*load_increment, 50 lz = clampMax[2]*load_increment, 51 // Normalized rotation axis 52 kx = clampMax[3], 53 ky = clampMax[4], 54 kz = clampMax[5], 55 // Rotation polynomial 56 c_0 = clampMax[6] * M_PI, 57 c_1 = clampMax[7] * M_PI, 58 cx = kx * x + ky * y + kz * z, 59 // Rotation magnitude 60 theta = (c_0 + c_1 * cx) * load_increment; 61 PetscScalar c = cos(theta), s = sin(theta); 62 63 u[0] = lx + s*(-kz*y + ky*z) + (1-c)*(-(ky*ky+kz*kz)*x + kx*ky*y + kx*kz*z); 64 u[1] = ly + s*(kz*x + -kx*z) + (1-c)*(kx*ky*x + -(kx*kx+kz*kz)*y + ky*kz*z); 65 u[2] = lz + s*(-ky*x + kx*y) + (1-c)*(kx*kz*x + ky*kz*y + -(kx*kx+ky*ky)*z); 66 PetscFunctionReturn(0); 67 }; 68