12 #ifndef OPM_ELASTICITY_IMPL_HPP
13 #define OPM_ELASTICITY_IMPL_HPP
18 namespace Elasticity {
20 template<
class Gr
idType>
21 template<
int components,
int funcdim>
23 const Dune::FieldVector<ctype,dim>& point,
24 const Dune::FieldMatrix<ctype,dim,dim>& Jinv)
28 int funcs = funcdim/dim;
30 Dune::FieldMatrix<
ctype,funcdim/dim,dim> dNdX;
31 for (
int i=0;i < funcs; i++)
32 Jinv.mv(basis[i].evaluateGradient(point),dNdX[i]);
33 static const int rows = 6+(dim-2)*12;
35 Dune::FieldMatrix<ctype,rows,funcdim/dim> temp;
39 #define INDEX(i,j) i+6*j
40 for (
int i=0; i < funcs; ++i) {
42 temp[
INDEX(0,0)][i] = dNdX[i][0];
43 temp[
INDEX(1,1)][i] = dNdX[i][1];
44 temp[
INDEX(2,2)][i] = dNdX[i][2];
47 temp[
INDEX(3,0)][i] = dNdX[i][1];
48 temp[
INDEX(3,1)][i] = dNdX[i][0];
50 temp[
INDEX(4,0)][i] = dNdX[i][2];
51 temp[
INDEX(4,2)][i] = dNdX[i][0];
52 temp[
INDEX(5,1)][i] = dNdX[i][2];
53 temp[
INDEX(5,2)][i] = dNdX[i][1];
55 }
else if (dim == 2) {
57 #define INDEX(i,j) i+3*j
58 for (
int i=0; i < funcs; ++i) {
60 temp[
INDEX(0,0)][i] = dNdX[i][0];
61 temp[
INDEX(1,1)][i] = dNdX[i][1];
64 temp[
INDEX(2,0)][i] = dNdX[i][1];
65 temp[
INDEX(2,1)][i] = dNdX[i][0];
69 for (
int j=0;
j<funcs*dim;++
j)
70 for (
size_t i=0;i<components;++i,++k)
71 B[i][
j] = temp[k%rows][k/rows];
74 template <
class Gr
idType>
75 template<
int comp,
int funcdim>
77 Dune::FieldMatrix<ctype,funcdim,funcdim>& A,
78 const Dune::FieldMatrix<ctype,comp,funcdim>& B,
79 const Dune::FieldMatrix<ctype,comp,comp>& C,
82 Dune::FieldMatrix<ctype,funcdim,comp> B2;
83 for (
int i=0;i<comp;++i)
84 for (
int j=0;
j<funcdim;++
j)
86 A = B.leftmultiplyany(C).leftmultiplyany(B2);
90 template<
class Gr
idType>
91 template<
int comp,
int funcdim>
93 const Dune::FieldVector<ctype,funcdim>& v,
94 const Dune::FieldVector<ctype,comp>&
eps0,
95 const Dune::FieldMatrix<ctype,comp,funcdim>& B,
96 const Dune::FieldMatrix<ctype,comp,comp>& C)
98 sigma = Dune::FMatrixHelp::mult(C,Dune::FMatrixHelp::mult(B,v)+eps0);
101 Dune::FieldVector<double,3>
waveSpeeds(
const Dune::FieldMatrix<double,6,6>& C,
double phi,
double theta,
double density)
104 Dune::FieldVector<double, 3> x;
105 x[0] = r*cos(theta)*cos(phi);
106 x[1] = r*sin(theta)*cos(phi);
109 Dune::FieldMatrix<double, 3, 6> D;
120 Dune::FieldMatrix<double, 6, 3> Dt;
121 for (
size_t i=0;i<6;++i)
122 for (
size_t j=0;
j<3;++
j)
125 Dune::FieldMatrix<double, 3, 6> T;
126 Dune::FieldMatrix<double, 3, 3> E;
127 Dune::FMatrixHelp::multMatrix(D, C, T);
128 Dune::FMatrixHelp::multMatrix(T, Dt, E);
129 Dune::FieldVector<double, 3> eigenvalues;
130 Dune::FMatrixHelp::eigenValues(E, eigenvalues);
131 std::sort(eigenvalues.begin(), eigenvalues.end(), std::greater<double>());
132 for (
size_t i=0;i<3;++i)
133 eigenvalues[i] = sqrt(eigenvalues[i]/density);
Definition: applier.hpp:18
Dune::FieldVector< double, 3 > waveSpeeds(const Dune::FieldMatrix< double, 6, 6 > &C, double phi, double theta, double density)
Compute the elastic wave velocities.
Definition: elasticity_impl.hpp:101
static const P1ShapeFunctionSet & instance()
Get the only instance of this class.
Definition: shapefunctions.hpp:193
Elasticity helper class.
Definition: elasticity.hpp:22
Dune::FieldVector< ctype, comp > eps0
Definition: elasticity_upscale_impl.hpp:445
GridType::LeafGridView::ctype ctype
A basic number.
Definition: elasticity.hpp:28
Classes for shape functions. Loosely based on code in dune-grid-howto.
Singleton handler for the set of LinearShapeFunction.
Definition: shapefunctions.hpp:180
int j
Definition: elasticity_upscale_impl.hpp:658