I am currently trying to solve the elastic wave equation using dG with symmetric interior penalties, and i seem to run into some problems regarding the discretization of the numerical flux. The flux-term i want to implement looks like
$$J=\int_\gamma \bigg([\phi]\{\sigma(u)\cdot n\}+\{\sigma(\phi)\cdot n\}[u]-\frac{\beta}{h}[\phi][u]\bigg) d\gamma$$ where $\phi$ would be a test field that is suitably nice, $u=[u_1,u_2]^T$ is the displacement field and $\sigma(u)$ the stress tensor. I have heard some people wanting to use the Voight notation, which seems reasonable but i cannot quite get it to work. I have tried writing the overall integramd above on the form
$$C^T\mathcal{T}\xi+\mathcal{T}^TC\xi-\frac{\beta}{h}C^TC\xi,$$ where $\xi=[\xi^+,\xi^-]^T$ are the degrees of freedom on the positive and negative element adjacent to $\gamma$ respectively, and
$$C=\begin{bmatrix}C^+ & 0 \\ 0 & -C^-\end{bmatrix}$$ is a matrix only containing the basis functions on each element, and $\mathcal{T}$ gives the Voight form of the stress tensor, i.e we have that
$$\mathcal{T}=\begin{bmatrix}NDB^+ & 0 \\ 0 & NDB^-\end{bmatrix}$$ with $D$ containing the Lamé parameters, $N=\begin{bmatrix}n_1 & 0 &n_2 \\ 0 & n_2 & n_1\end{bmatrix}$ a matrix containing the outward unit normal from the positive element and $B$ a matrix containing the derivatives of the basis functions so that
$$B^\pm\xi^\pm=\varepsilon_v^\pm=(\varepsilon^\pm_1,\varepsilon^\pm_2,2\varepsilon^\pm_{12})^T.$$
With this formulation i have also run into the problem of how the components of the total matrix in the integrand should be sent out to the global flux matrix. What would be a clear way of doing this? I feel like there arises some difficulties if you e.g. consider the degrees of freedom be ordered at each element as $\xi^+=[\xi^+_x,\xi^+_y]^T$ which i want to have in order to be compatible with another code i want to use this with later on.
So to summarize the questions:
- Is the Voight method best? Or do you use any other way of evaluating the integrand?
- How do you map it to the global flux matrix in an efficient manner
Thanks in advance!