diff --git a/src/rodassembler.cc b/src/rodassembler.cc index 2827907480cf10a382ab0dd67e6fb9fed491f928..bb09d1e1f9a2cddb8a0683fe336469722f4730b5 100644 --- a/src/rodassembler.cc +++ b/src/rodassembler.cc @@ -158,8 +158,8 @@ getStrain(const std::vector<Configuration>& localSolution, template <class GridType> -void Dune::RodAssembler<GridType>:: -getNeighborsPerVertex(MatrixIndexSet& nb) const +void RodAssembler<GridType>:: +getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const { const int gridDim = GridType::dimension; const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel()); @@ -193,10 +193,12 @@ getNeighborsPerVertex(MatrixIndexSet& nb) const template <class GridType> -void Dune::RodAssembler<GridType>:: +void RodAssembler<GridType>:: assembleMatrix(const std::vector<Configuration>& sol, - BCRSMatrix<MatrixBlock>& matrix) const + Dune::BCRSMatrix<MatrixBlock>& matrix) const { + using namespace Dune; + const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel()); MatrixIndexSet neighborsPerVertex; @@ -246,12 +248,14 @@ assembleMatrix(const std::vector<Configuration>& sol, template <class GridType> template <class MatrixType> -void Dune::RodAssembler<GridType>:: +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; + /* ndof is the number of vectors of the element */ int ndof = matSize; @@ -528,10 +532,12 @@ getLocalMatrix( EntityPointer &entity, } template <class GridType> -void Dune::RodAssembler<GridType>:: +void RodAssembler<GridType>:: assembleMatrixFD(const std::vector<Configuration>& sol, - BCRSMatrix<MatrixBlock>& matrix) const + Dune::BCRSMatrix<MatrixBlock>& matrix) const { + using namespace Dune; + double eps = 1e-2; typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> >::row_type::iterator ColumnIterator; @@ -730,23 +736,16 @@ assembleMatrixFD(const std::vector<Configuration>& sol, } template <class GridType> -void Dune::RodAssembler<GridType>:: +void RodAssembler<GridType>:: strainDerivative(const std::vector<Configuration>& localSolution, double pos, - FieldVector<double,1> shapeGrad[2], - FieldVector<double,1> shapeFunction[2], - Dune::array<FieldMatrix<double,2,6>, 6>& derivatives) const + Dune::FieldVector<double,1> shapeGrad[2], + Dune::FieldVector<double,1> shapeFunction[2], + Dune::array<Dune::FieldMatrix<double,2,6>, 6>& derivatives) const { - assert(localSolution.size()==2); - -// FieldVector<double,1> shapeGrad[2]; -// shapeGrad[0] = -1; -// shapeGrad[1] = 1; - -// FieldVector<double,1> shapeFunction[2]; -// shapeFunction[0] = 1-pos; -// shapeFunction[1] = pos; + using namespace Dune; + assert(localSolution.size()==2); FieldVector<double,3> r_s; for (int i=0; i<3; i++) @@ -814,12 +813,14 @@ strainDerivative(const std::vector<Configuration>& localSolution, template <class GridType> -void Dune::RodAssembler<GridType>:: +void RodAssembler<GridType>:: strainHessian(const std::vector<Configuration>& localSolution, double pos, - Dune::array<Matrix<FieldMatrix<double,6,6> >, 3>& translationDer, - Dune::array<Matrix<FieldMatrix<double,3,3> >, 3>& rotationDer) const + Dune::array<Dune::Matrix<Dune::FieldMatrix<double,6,6> >, 3>& translationDer, + Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3>& rotationDer) const { + using namespace Dune; + assert(localSolution.size()==2); FieldVector<double,1> shapeGrad[2]; @@ -942,13 +943,15 @@ strainHessian(const std::vector<Configuration>& localSolution, } template <class GridType> -void Dune::RodAssembler<GridType>:: +void RodAssembler<GridType>:: rotationStrainHessian(const std::vector<Configuration>& x, double pos, Dune::FieldVector<double,1> shapeGrad[2], Dune::FieldVector<double,1> shapeFunction[2], Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3>& rotationDer) const { + using namespace Dune; + assert(x.size()==2); double eps = 1e-3; @@ -1003,11 +1006,13 @@ rotationStrainHessian(const std::vector<Configuration>& x, template <class GridType> -void Dune::RodAssembler<GridType>:: +void RodAssembler<GridType>:: assembleGradient(const std::vector<Configuration>& sol, - BlockVector<FieldVector<double, blocksize> >& grad) const + Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const { - const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel()); + using namespace Dune; + + const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel()); const int maxlevel = grid_->maxLevel(); if (sol.size()!=grid_->size(maxlevel, gridDim)) @@ -1180,9 +1185,11 @@ assembleGradient(const std::vector<Configuration>& sol, template <class GridType> -double Dune::RodAssembler<GridType>:: +double RodAssembler<GridType>:: computeEnergy(const std::vector<Configuration>& sol) const { + using namespace Dune; + double energy = 0; const typename GridType::Traits::LeafIndexSet& indexSet = grid_->leafIndexSet(); @@ -1268,10 +1275,12 @@ computeEnergy(const std::vector<Configuration>& sol) const template <class GridType> -void Dune::RodAssembler<GridType>:: +void RodAssembler<GridType>:: getStrain(const std::vector<Configuration>& sol, - BlockVector<FieldVector<double, blocksize> >& strain) const + Dune::BlockVector<Dune::FieldVector<double, blocksize> >& strain) const { + using namespace Dune; + const typename GridType::Traits::LeafIndexSet& indexSet = grid_->leafIndexSet(); if (sol.size()!=indexSet.size(gridDim)) @@ -1330,10 +1339,12 @@ getStrain(const std::vector<Configuration>& sol, } template <class GridType> -Dune::FieldVector<double, 6> Dune::RodAssembler<GridType>::getStrain(const std::vector<Configuration>& sol, +Dune::FieldVector<double, 6> RodAssembler<GridType>::getStrain(const std::vector<Configuration>& sol, const EntityPointer& element, double pos) const { + using namespace Dune; + if (!element->isLeaf()) DUNE_THROW(Dune::NotImplemented, "Only for leaf elements"); @@ -1416,11 +1427,13 @@ Dune::FieldVector<double, 6> Dune::RodAssembler<GridType>::getStrain(const std:: } template <class GridType> -Dune::FieldVector<double,3> Dune::RodAssembler<GridType>:: +Dune::FieldVector<double,3> RodAssembler<GridType>:: getResultantForce(const BoundaryPatch<GridType>& boundary, const std::vector<Configuration>& sol, - FieldVector<double,3>& canonicalTorque) const + Dune::FieldVector<double,3>& canonicalTorque) const { + using namespace Dune; + const typename GridType::Traits::LeafIndexSet& indexSet = grid_->leafIndexSet(); if (sol.size()!=indexSet.size(gridDim)) diff --git a/src/rodassembler.hh b/src/rodassembler.hh index fe98940a3b0d2a115299f7cac133ea91806c2fee..d35dbc516eae56579621a4708f6e7c77ea591477 100644 --- a/src/rodassembler.hh +++ b/src/rodassembler.hh @@ -117,8 +117,6 @@ public: }; -namespace Dune -{ /** \brief The FEM operator for an extensible, shearable rod */ @@ -139,7 +137,7 @@ namespace Dune enum { blocksize = 6 }; //! - typedef FieldMatrix<double, blocksize, blocksize> MatrixBlock; + typedef Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock; /** \todo public only for debugging! */ public: @@ -221,60 +219,60 @@ namespace Dune A_[1] = G * A; A_[2] = E * A; - printf("%g %g %g %g %g %g\n", K_[0], K_[1], K_[2], A_[0], A_[1], A_[2]); + //printf("%g %g %g %g %g %g\n", K_[0], K_[1], K_[2], A_[0], A_[1], A_[2]); //exit(0); } /** \brief Assemble the tangent stiffness matrix */ void assembleMatrix(const std::vector<Configuration>& sol, - BCRSMatrix<MatrixBlock>& matrix) const; + Dune::BCRSMatrix<MatrixBlock>& matrix) const; /** \brief Assemble the tangent stiffness matrix using a finite difference approximation */ void assembleMatrixFD(const std::vector<Configuration>& sol, - BCRSMatrix<MatrixBlock>& matrix) const; + Dune::BCRSMatrix<MatrixBlock>& matrix) const; void strainDerivative(const std::vector<Configuration>& localSolution, double pos, - FieldVector<double,1> shapeGrad[2], - FieldVector<double,1> shapeFunction[2], - Dune::array<FieldMatrix<double,2,6>, 6>& derivatives) const; + Dune::FieldVector<double,1> shapeGrad[2], + Dune::FieldVector<double,1> shapeFunction[2], + Dune::array<Dune::FieldMatrix<double,2,6>, 6>& derivatives) const; void rotationStrainHessian(const std::vector<Configuration>& x, double pos, - FieldVector<double,1> shapeGrad[2], - FieldVector<double,1> shapeFunction[2], + Dune::FieldVector<double,1> shapeGrad[2], + Dune::FieldVector<double,1> shapeFunction[2], Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3>& rotationDer) const; void strainHessian(const std::vector<Configuration>& localSolution, double pos, - Dune::array<Matrix<FieldMatrix<double,6,6> >, 3>& translationDer, - Dune::array<Matrix<FieldMatrix<double,3,3> >, 3>& rotationDer) const; + Dune::array<Dune::Matrix<Dune::FieldMatrix<double,6,6> >, 3>& translationDer, + Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3>& rotationDer) const; void assembleGradient(const std::vector<Configuration>& sol, - BlockVector<FieldVector<double, blocksize> >& grad) const; + Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const; /** \brief Compute the energy of a deformation state */ double computeEnergy(const std::vector<Configuration>& sol) const; - void getNeighborsPerVertex(MatrixIndexSet& nb) const; + void getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const; void getStrain(const std::vector<Configuration>& sol, - BlockVector<FieldVector<double, blocksize> >& strain) const; + Dune::BlockVector<Dune::FieldVector<double, blocksize> >& strain) const; /** \brief Get the strain at a particular point of the grid */ - FieldVector<double, 6> getStrain(const std::vector<Configuration>& sol, - const EntityPointer& element, - double pos) const; + Dune::FieldVector<double, 6> getStrain(const std::vector<Configuration>& sol, + const EntityPointer& element, + double pos) const; /** \brief Return resultant force across boundary in canonical coordinates \note Linear run-time in the size of the grid */ - FieldVector<double,3> getResultantForce(const BoundaryPatch<GridType>& boundary, - const std::vector<Configuration>& sol, - FieldVector<double,3>& canonicalTorque) const; + Dune::FieldVector<double,3> getResultantForce(const BoundaryPatch<GridType>& boundary, + const std::vector<Configuration>& sol, + Dune::FieldVector<double,3>& canonicalTorque) const; protected: @@ -287,9 +285,9 @@ namespace Dune const int matSize, MatrixType& mat) const; template <class T> - static FieldVector<T,3> darboux(const Quaternion<T>& q, const FieldVector<T,4>& q_s) + static Dune::FieldVector<T,3> darboux(const Quaternion<T>& q, const Dune::FieldVector<T,4>& q_s) { - FieldVector<double,3> u; // The Darboux vector + Dune::FieldVector<double,3> u; // The Darboux vector u[0] = 2 * (q.B(0) * q_s); u[1] = 2 * (q.B(1) * q_s); @@ -299,8 +297,6 @@ namespace Dune } }; // end class - -} // end namespace #include "rodassembler.cc" diff --git a/src/rodsolver.cc b/src/rodsolver.cc index b3d454118af09cd7e274bca82e2a612a286d708e..869f4549cfc0e2320d71dda61b9db8def7fd6058 100644 --- a/src/rodsolver.cc +++ b/src/rodsolver.cc @@ -47,7 +47,7 @@ setTrustRegionObstacles(double trustRegionRadius, template <class GridType> void RodSolver<GridType>::setup(const GridType& grid, - const Dune::RodAssembler<GridType>* rodAssembler, + const RodAssembler<GridType>* rodAssembler, const SolutionType& x, int maxTrustRegionSteps, double initialTrustRegionRadius, diff --git a/src/rodsolver.hh b/src/rodsolver.hh index a2a11dd9ff485c8dd4361c1ec9b829c4097c1f73..4b8dcc496b87fba2daad1ac3754170127c24672e 100644 --- a/src/rodsolver.hh +++ b/src/rodsolver.hh @@ -33,7 +33,7 @@ public: {} void setup(const GridType& grid, - const Dune::RodAssembler<GridType>* rodAssembler, + const RodAssembler<GridType>* rodAssembler, const SolutionType& x, int maxTrustRegionSteps, double initialTrustRegionRadius, @@ -92,7 +92,7 @@ protected: MatrixType* hessianMatrix_; /** \brief The assembler for the material law */ - const Dune::RodAssembler<GridType>* rodAssembler_; + const RodAssembler<GridType>* rodAssembler_; /** \brief The multigrid solver */ Dune::IterativeSolver<MatrixType, CorrectionType>* mmgSolver_;