Go to the documentation of this file.
12#ifndef OPM_ELASTICITY_IMPL_HPP
13#define OPM_ELASTICITY_IMPL_HPP
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);
101Dune::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);
Elasticity helper class. Definition: elasticity.hpp:22
GridType::LeafGridView::ctype ctype A basic number. Definition: elasticity.hpp:28
Singleton handler for the set of LinearShapeFunction. Definition: shapefunctions.hpp:181
static const P1ShapeFunctionSet & instance() Get the only instance of this class. Definition: shapefunctions.hpp:193
int j Definition: elasticity_upscale_impl.hpp:658
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
Dune::FieldVector< ctype, comp > eps0 Definition: elasticity_upscale_impl.hpp:446
Definition: ImplicitAssembly.hpp:43
Classes for shape functions. Loosely based on code in dune-grid-howto.
|