diff --git a/src/rodassembler.cc b/src/rodassembler.cc
index 969fa4f787958615752c5aa12406669ff2d2a8da..2983afcf91d02aba121ea0176fa4c80491699e6b 100644
--- a/src/rodassembler.cc
+++ b/src/rodassembler.cc
@@ -561,8 +561,9 @@ assembleMatrix(const std::vector<Configuration>& sol,
localSolution[i] = sol[indexSet.template subIndex<gridDim>(*it,i)];
// setup matrix
- getLocalMatrix( it, localSolution, sol, numOfBaseFct, mat);
-
+ //getLocalMatrix( it, localSolution, sol, numOfBaseFct, mat);
+ DUNE_THROW(NotImplemented, "getLocalMatrix");
+
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
@@ -580,292 +581,6 @@ assembleMatrix(const std::vector<Configuration>& sol,
}
-template <class GridType>
-template <class MatrixType>
-void RodAssembler<GridType>::
-getLocalMatrix( EntityPointer &entity,
- const std::vector<Configuration>& localSolution,
- const std::vector<Configuration>& globalSolution,
- const int matSize, MatrixType& localMat) const
-{
- using namespace Dune;
-#if 0
-
- /* ndof is the number of vectors of the element */
- int ndof = matSize;
-
- for (int i=0; i<matSize; i++)
- for (int j=0; j<matSize; j++)
- localMat[i][j] = 0;
-
- const LagrangeShapeFunctionSet<double, double, gridDim> & baseSet
- = Dune::LagrangeShapeFunctions<double, double, gridDim>::general(entity->type(), elementOrder);
-
- // Get quadrature rule
- int polOrd = 2;
- const QuadratureRule<double, gridDim>& quad = QuadratureRules<double, gridDim>::rule(entity->type(), polOrd);
-
- /* Loop over all integration points */
- for (size_t ip=0; ip<quad.size(); ip++) {
-
- // Local position of the quadrature point
- const FieldVector<double,gridDim>& quadPos = quad[ip].position();
-
- // calc Jacobian inverse before integration element is evaluated
- const FieldMatrix<double,gridDim,gridDim>& inv = entity->geometry().jacobianInverseTransposed(quadPos);
- const double integrationElement = entity->geometry().integrationElement(quadPos);
-
- /* Compute the weight of the current integration point */
- double weight = quad[ip].weight() * integrationElement;
-
- /**********************************************/
- /* compute gradients of the shape functions */
- /**********************************************/
- FieldVector<double,gridDim> shapeGrad[2];
- FieldVector<double,gridDim> tmp;
-
- for (int dof=0; dof<ndof; dof++) {
-
- for (int i=0; i<gridDim; i++)
- tmp[i] = baseSet[dof].evaluateDerivative(0,i,quadPos);
-
- // multiply with jacobian inverse
- shapeGrad[dof] = 0;
- inv.umv(tmp,shapeGrad[dof]);
-
- }
-
-
- FieldVector<double,1> shapeFunction[matSize];
- for(int i=0; i<matSize; i++)
- shapeFunction[i] = baseSet[i].evaluateFunction(0,quadPos);
-
- // //////////////////////////////////
- // Interpolate
- // //////////////////////////////////
-
- FieldVector<double,3> r_s;
- r_s[0] = localSolution[0].r[0]*shapeGrad[0] + localSolution[1].r[0]*shapeGrad[1];
- r_s[1] = localSolution[0].r[1]*shapeGrad[0] + localSolution[1].r[1]*shapeGrad[1];
- r_s[2] = localSolution[0].r[2]*shapeGrad[0] + localSolution[1].r[2]*shapeGrad[1];
-
- // Interpolate current rotation at this quadrature point
- Quaternion<double> q = Quaternion<double>::interpolate(localSolution[0].q, localSolution[1].q,quadPos[0]);
-
- // Contains \partial q / \partial v^i_j at v = 0
- array<Quaternion<double>,3> dq_dvj;
- Quaternion<double> dq_dvij_ds[2][3];
- for (int i=0; i<2; i++)
- for (int j=0; j<3; j++) {
-
- for (int m=0; m<3; m++) {
- dq_dvj[j][m] = (j==m) * 0.5;
- dq_dvij_ds[i][j][m] = (j==m) * 0.5 * shapeGrad[i];
- }
-
- dq_dvj[j][3] = 0;
- dq_dvij_ds[i][j][3] = 0;
- }
-
- Quaternion<double> dq_dvj_dvl[3][3];
- Quaternion<double> dq_dvij_dvkl_ds[2][3][2][3];
- for (int i=0; i<2; i++) {
-
- for (int j=0; j<3; j++) {
-
- for (int k=0; k<2; k++) {
-
- for (int l=0; l<3; l++) {
-
- for (int m=0; m<3; m++) {
- dq_dvj_dvl[j][l][m] = 0;
- dq_dvij_dvkl_ds[i][j][k][l][m] = 0;
- }
-
- dq_dvj_dvl[j][l][3] = -0.25 * (j==l);
- dq_dvij_dvkl_ds[i][j][k][l][3] = -0.25 * (j==l) * shapeGrad[i] * shapeGrad[k];
-
- }
-
- }
-
- }
-
- }
-
- // Contains \parder d \parder v^i_j
- array<array<FieldVector<double,3>, 3>, 3> dd_dvj;
- q.getFirstDerivativesOfDirectors(dd_dvj);
-
- // Contains \parder {dm}{v^i_j}{v^k_l}
- FieldVector<double,3> dd_dvij_dvkl[3][3][3];
-
- for (int j=0; j<3; j++) {
-
- for (int l=0; l<3; l++) {
-
- FieldMatrix<double,4,4> A;
- for (int a=0; a<4; a++)
- for (int b=0; b<4; b++)
- A[a][b] = (q.mult(dq_dvj[l]))[a] * (q.mult(dq_dvj[j]))[b]
- + q[a] * q.mult(dq_dvj_dvl[l][j])[b];
-
- // d1
- dd_dvij_dvkl[0][j][l][0] = A[0][0] - A[1][1] - A[2][2] + A[3][3];
- dd_dvij_dvkl[0][j][l][1] = A[1][0] + A[0][1] + A[3][2] + A[2][3];
- dd_dvij_dvkl[0][j][l][2] = A[2][0] + A[0][2] - A[3][1] - A[1][3];
-
- // d2
- dd_dvij_dvkl[1][j][l][0] = A[1][0] + A[0][1] - A[3][2] - A[2][3];
- dd_dvij_dvkl[1][j][l][1] = -A[0][0] + A[1][1] - A[2][2] + A[3][3];
- dd_dvij_dvkl[1][j][l][2] = A[2][1] + A[1][2] + A[3][0] + A[0][3];
-
- // d3
- dd_dvij_dvkl[2][j][l][0] = A[2][0] + A[0][2] + A[3][1] + A[1][3];
- dd_dvij_dvkl[2][j][l][1] = A[2][1] + A[1][2] - A[3][0] - A[0][3];
- dd_dvij_dvkl[2][j][l][2] = -A[0][0] - A[1][1] + A[2][2] + A[3][3];
-
-
- dd_dvij_dvkl[0][j][l] *= 2;
- dd_dvij_dvkl[1][j][l] *= 2;
- dd_dvij_dvkl[2][j][l] *= 2;
-
- }
-
- }
-
- // Get the derivative of the rotation at the quadrature point by interpolating in $H$
- Quaternion<double> q_s = Quaternion<double>::interpolateDerivative(localSolution[0].q, localSolution[1].q,
- quadPos, 1/shapeGrad[1]);
-
- // The current strain
- FieldVector<double,blocksize> strain = getStrain(globalSolution, entity, quadPos);
-
- // The reference strain
- FieldVector<double,blocksize> referenceStrain = getStrain(referenceConfiguration_, entity, quadPos);
-
- // Contains \partial u (canonical) / \partial v^i_j at v = 0
- FieldVector<double,3> du_dvij[2][3];
- FieldVector<double,3> du_dvij_dvkl[2][3][2][3];
-
- for (int i=0; i<2; i++) {
-
- for (int j=0; j<3; j++) {
-
- for (int m=0; m<3; m++) {
- du_dvij[i][j][m] = (q.mult(dq_dvj[j])).B(m)*q_s;
- du_dvij[i][j][m] *= 2 * shapeFunction[i];
- Quaternion<double> tmp = dq_dvj[j];
- tmp *= shapeFunction[i];
- du_dvij[i][j][m] += 2 * ( q.B(m)*(q.mult(dq_dvij_ds[i][j]) + q_s.mult(tmp)));
- }
-
-#if 1
- for (int k=0; k<2; k++) {
-
- for (int l=0; l<3; l++) {
-
- Quaternion<double> tmp_ij = dq_dvj[j];
- Quaternion<double> tmp_kl = dq_dvj[l];
-
- tmp_ij *= shapeFunction[i];
- tmp_kl *= shapeFunction[k];
-
- Quaternion<double> tmp_ijkl = dq_dvj_dvl[j][l];
- tmp_ijkl *= shapeFunction[i] * shapeFunction[k];
-
- for (int m=0; m<3; m++) {
-
- du_dvij_dvkl[i][j][k][l][m] = 2 * ( (q.mult(tmp_ijkl)).B(m) * q_s)
- + 2 * ( (q.mult(tmp_ij)).B(m) * (q.mult(dq_dvij_ds[k][l]) + q_s.mult(tmp_kl)))
- + 2 * ( (q.mult(tmp_kl)).B(m) * (q.mult(dq_dvij_ds[i][j]) + q_s.mult(tmp_ij)))
- + 2 * ( q.B(m) * (q.mult(dq_dvij_dvkl_ds[i][j][k][l]) + q_s.mult(tmp_ijkl)));
-
- }
-
- }
-
- }
-#else
-#warning Using fd approximation of rotation strain Hessian.
- Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3> rotationHessian;
- rotationStrainHessian(localSolution, quadPos[0], shapeGrad, shapeFunction, rotationHessian);
- for (int k=0; k<2; k++)
- for (int l=0; l<3; l++)
- for (int m=0; m<3; m++)
- du_dvij_dvkl[i][j][k][l][m] = rotationHessian[m][i][k][j][l];
-#endif
-
- }
-
- }
-
- // ///////////////////////////////////
- // Sum it all up
- // ///////////////////////////////////
- array<FieldMatrix<double,2,6>, 6> strainDerivatives;
- strainDerivative(localSolution, quadPos[0], shapeGrad, shapeFunction, strainDerivatives);
-
- for (int i=0; i<matSize; i++) {
-
- for (int k=0; k<matSize; k++) {
-
- for (int j=0; j<3; j++) {
-
- for (int l=0; l<3; l++) {
-
- for (int m=0; m<3; m++) {
-
- // ////////////////////////////////////////////
- // The translational part
- // ////////////////////////////////////////////
-
- // \partial W^2 \partial r^i_j \partial r^k_l
- localMat[i][k][j][l] += weight
- * A_[m] * shapeGrad[i] * q.director(m)[j] * shapeGrad[k] * q.director(m)[l];
-
- // \partial W^2 \partial v^i_j \partial r^k_l
- localMat[i][k][j+3][l] += weight
- * ( A_[m] * strainDerivatives[m][k][l] * strainDerivatives[m][i][j+3]
- + A_[m] * (strain[m] - referenceStrain[m]) * shapeGrad[k] * dd_dvj[m][j][l]*shapeFunction[i]);
-
- // \partial W^2 \partial r^i_j \partial v^k_l
- localMat[i][k][j][l+3] += weight
- * ( A_[m] * strainDerivatives[m][i][j] * strainDerivatives[m][k][l+3]
- + A_[m] * (strain[m] - referenceStrain[m]) * shapeGrad[i] * dd_dvj[m][l][j]*shapeFunction[k]);
-
- // \partial W^2 \partial v^i_j \partial v^k_l
- localMat[i][k][j+3][l+3] += weight
- * ( A_[m] * (r_s * dd_dvj[m][l]*shapeFunction[k]) * (r_s * dd_dvj[m][j]*shapeFunction[i])
- + A_[m] * (strain[m] - referenceStrain[m]) * (r_s * dd_dvij_dvkl[m][j][l]) * shapeFunction[i] * shapeFunction[k]);
-
- // ////////////////////////////////////////////
- // The rotational part
- // ////////////////////////////////////////////
-
- // \partial W^2 \partial v^i_j \partial v^k_l
- // All other derivatives are zero
-
- double sum = du_dvij[k][l][m] * du_dvij[i][j][m];
-
- sum += (strain[m+3] - referenceStrain[m+3]) * du_dvij_dvkl[i][j][k][l][m];
-
- localMat[i][k][j+3][l+3] += weight *K_[m] * sum;
-
- }
-
- }
-
- }
-
- }
-
- }
-
- }
-#endif
-}
-
template <class GridType>
void RodAssembler<GridType>::
assembleMatrixFD(const std::vector<Configuration>& sol,
diff --git a/src/rodassembler.hh b/src/rodassembler.hh
index 51d1cb4b67311fc4f1565d1ea68ee94b9f6bac48..3b4dda9bf6f6a691550a2c0489881d18acc2e1e7 100644
--- a/src/rodassembler.hh
+++ b/src/rodassembler.hh
@@ -285,14 +285,6 @@ public:
protected:
- /** \brief Compute the element tangent stiffness matrix
- \todo Handing over both the local and the global solution is pretty stupid. */
- template <class MatrixType>
- void getLocalMatrix( EntityPointer &entity,
- const std::vector<Configuration>& localSolution,
- const std::vector<Configuration>& globalSolution,
- const int matSize, MatrixType& mat) const;
-
template <class T>
static Dune::FieldVector<T,3> darboux(const Quaternion<T>& q, const Dune::FieldVector<T,4>& q_s)
{