dune-localfunctions
2.2.0
|
00001 // -*- tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=2 sw=2 sts=2: 00003 #ifndef DUNE_Q2_LOCALBASIS_HH 00004 #define DUNE_Q2_LOCALBASIS_HH 00005 00006 #include <dune/common/fmatrix.hh> 00007 00008 #include <dune/localfunctions/common/localbasis.hh> 00009 00010 namespace Dune 00011 { 00023 template<class D, class R, int dim> 00024 class Q2LocalBasis 00025 { 00026 public: 00027 typedef LocalBasisTraits<D,dim,Dune::FieldVector<D,dim>,R,1,Dune::FieldVector<R,1>, 00028 Dune::FieldMatrix<R,1,dim> > Traits; 00029 00031 unsigned int size () const 00032 { 00033 int size = 1; 00034 for (int i=0; i<dim; i++) 00035 size *= 3; 00036 return size; 00037 } 00038 00040 inline void evaluateFunction (const typename Traits::DomainType& in, 00041 std::vector<typename Traits::RangeType>& out) const 00042 { 00043 out.resize(size()); 00044 00045 // Evaluate the Lagrange functions 00046 array<array<R,3>, dim> X; 00047 00048 for (size_t i=0; i<dim; i++) { 00049 X[i][0] = R(2)*in[i]*in[i] - R(3)*in[i]+R(1); 00050 X[i][1] = -R(4)*in[i]*in[i] + R(4)*in[i]; 00051 X[i][2] = R(2)*in[i]*in[i] - in[i]; 00052 } 00053 00054 00055 // Compute function values: they are products of 1d Lagrange function values 00056 for (size_t i=0; i<out.size(); i++) { 00057 00058 out[i] = 1; 00059 00060 // Construct the i-th Lagrange point 00061 size_t ternary = i; 00062 for (int j=0; j<dim; j++) { 00063 00064 int digit = ternary%3; 00065 ternary /= 3; 00066 00067 // Multiply the 1d Lagrange shape functions together 00068 out[i] *= X[j][digit]; 00069 00070 } 00071 00072 } 00073 00074 } 00075 00077 inline void 00078 evaluateJacobian (const typename Traits::DomainType& in, // position 00079 std::vector<typename Traits::JacobianType>& out) const // return value 00080 { 00081 out.resize(size()); 00082 00083 // Evaluate the 1d Lagrange functions and their derivatives 00084 array<array<R,3>, dim> X, DX; 00085 00086 for (size_t i=0; i<dim; i++) { 00087 X[i][0] = R(2)*in[i]*in[i] - R(3)*in[i]+R(1); 00088 X[i][1] = -R(4)*in[i]*in[i] + R(4)*in[i]; 00089 X[i][2] = R(2)*in[i]*in[i] - in[i]; 00090 00091 DX[i][0] = R(4)*in[i] - R(3); 00092 DX[i][1] = -R(8)*in[i] + R(4); 00093 DX[i][2] = R(4)*in[i] - R(1); 00094 } 00095 00096 00097 // Compute the derivatives by deriving the products of 1d Lagrange functions 00098 for (size_t i=0; i<out.size(); i++) { 00099 00100 // Computing the j-th partial derivative 00101 for (int j=0; j<dim; j++) { 00102 00103 out[i][0][j] = 1; 00104 00105 // Loop over the 'dim' terms in the product rule 00106 size_t ternary = i; 00107 for (int k=0; k<dim; k++) { 00108 00109 int digit = ternary%3; 00110 ternary /= 3; 00111 00112 out[i][0][j] *= (k==j) ? DX[k][digit] : X[k][digit]; 00113 00114 } 00115 00116 } 00117 00118 } 00119 00120 00121 00122 } 00123 00125 unsigned int order () const 00126 { 00127 return 2; 00128 } 00129 }; 00130 } 00131 #endif