The Method of manufactured is commonly used for verification of computational science codes. I want to use the method for verification of Navier-Cauchy (elasticity) equations with periodic and Bloch-periodic boundary conditions. In this case, the coefficients of the PDE (material propertise) are periodic, and the form of the boundary value problem is
$$\nabla\cdot \sigma + \mathbf{f} = -\omega^2 \rho \mathbf{u}\quad \forall \mathbf{x} \in \Omega\\ \mathbf{u}(\mathbf{x} + \mathbf{a}) = \mathbf{u}(\mathbf{x})e^{\mathbf{k}\cdot\mathbf{a}}\quad \forall \mathbf{x} \in \Gamma_u\\ \mathbf{t}(\mathbf{x} + \mathbf{a}) = -\mathbf{t}(\mathbf{x})e^{\mathbf{k}\cdot\mathbf{a}}\quad \forall \mathbf{x} \in \Gamma_t $$ with $$\begin{align} &\sigma = \underline{C} \epsilon\\ &\epsilon = \frac{1}{2}\left[\nabla\mathbf{u} + \nabla\mathbf{u}^T\right] \enspace , \end{align}$$ and $\underline{C}$ presents the desired periodicity, i.e., $\underline{C}(\mathbf{x} + \mathbf{a}) = \underline{C}(\mathbf{x})$, being $\mathbf{a}$ the vector periodicity of the material.
When $\omega=0$ the Bloch-periodic conditions turn into common periodic conditions and $\mathbf{u}(\mathbf{x} + \mathbf{a}) = \mathbf{u}(\mathbf{x})$ (similar for the tractions $\mathbf{t}(\mathbf{x})$). In that case one can propose a solution that presents the same periodicity of the material, e.g
$$\mathbf{u} = \left[e^{\sin \left(2\,x\right)\,\sin y}-\cos \left(3\,x\right)\,\sin ^2 \left(2\,y\right)\right] (u_1, u_2, 0)$$ and the periodicity for the material being $$\underline{C} = \left[\sin ^2x\,\sin \cos y^3+\sin \cos x^2\,\sin ^4y\right] \underline{C}_0 \enspace ,$$
then, plug these terms in the equations and obtain the expression for $\mathbf{f}$ (this can be done symbolically using Maxima or SymPy and then generate Fortran, C++ or Python code).
In the case of Bloch-periodicity it is a little bit more convoluted, since the functions that solve the equation are not periodic. They are the product of two periodic functions, and just periodic in the particular case of commensurate wavelengths. Thus, in this case the solution would depend on the wavevector $\mathbf{k}$.
Question: It is customary to assume that the material properties are composites where the properties change abruptly, e.g., a bilayer material, fiber reinforced composites. In these cases the constitutive tensor is not a differentiable function. If the code is verified with smooth functions what can be said about the case of non-differentiable material properties? Is the test still valid?