Method of Manufactured Solutions for non-differentiable coefficients

Question

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?

Answer 1

Let me see if I can turn my comment into a fully-fledged answer. Suppose you are solving the Laplacian with a variable $k$ which has a step function across the domain. Suppose you have homogeneous boundary conditions so that we don't have to deal with that in my answer. Now your weak form will look something like:

$$ \int_\Omega k\nabla u \cdot \nabla v \; dx = \int_\Omega f v \;dx \;\;\forall v $$

Suppose you want to have $u=\tanh(x)\sin(x)$, in the strong form, you can't differentiate $k$, because it's discontinuous, but you can compute

$$ \int_\Omega k\nabla (\tanh(x)\sin(x)) \cdot \nabla v \; dx = \int_\Omega f v \;dx \;\;\forall v $$

Now, given that, instead of trying to compute $f$ from this, you can just implement the above as your right-hand side. I.e., solve

$$ \int_\Omega k\nabla u_h \cdot \nabla v_h \; dx = \int_\Omega k\nabla (\tanh(x)\sin(x)) \cdot \nabla v_h \; dx \;\;\forall v_h $$

using your favorite $k$. As long as your integration makes sense, and as long as you can take that single derivative of your chosen solution, you should be good to go.