diff --git a/rodobstacle.cc b/rodobstacle.cc index c07ba497d8b0eed24f40208ad6f524abbd2e977b..c6938d88e00414207a7d2943d6f0ee6dd8007c28 100644 --- a/rodobstacle.cc +++ b/rodobstacle.cc @@ -152,11 +152,12 @@ int main (int argc, char *argv[]) try // ////////////////////////// for (int i=0; i<x.size(); i++) { - x[i].r[1] = i;//double(i)/(x.size()-1); x[i].r[0] = 0; + x[i].r[1] = i;//double(i)/(x.size()-1); x[i].q = Rotation<2,double>::identity(); } + x.back().r[1] += 1; // ///////////////////////////////////////////////////////////////////// // Refinement Loop @@ -180,7 +181,7 @@ int main (int argc, char *argv[]) try MatrixType hessianMatrix; - PlanarRodAssembler<GridType,4> rodAssembler(grid); + PlanarRodAssembler<GridType> rodAssembler(grid); rodAssembler.setParameters(1, 350000, 350000); diff --git a/src/rodassembler.cc b/src/rodassembler.cc index 5e439209e3d41f85f1b94dc15ac3bed6b375204d..db05fafb3bdcfdd0491047a78825ec088f4a3b45 100644 --- a/src/rodassembler.cc +++ b/src/rodassembler.cc @@ -214,8 +214,8 @@ getResultantForce(const BoundaryPatchBase<PatchGridView>& boundary, } -template <class GridType, int polOrd> -void PlanarRodAssembler<GridType, polOrd>::getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const +template <class GridType> +void PlanarRodAssembler<GridType>::getNeighborsPerVertex(Dune::MatrixIndexSet& nb) const { const int gridDim = GridType::dimension; const typename GridType::Traits::LevelIndexSet& indexSet = grid_->levelIndexSet(grid_->maxLevel()); @@ -247,8 +247,8 @@ void PlanarRodAssembler<GridType, polOrd>::getNeighborsPerVertex(Dune::MatrixInd } -template <class GridType, int polOrd> -void PlanarRodAssembler<GridType, polOrd>:: +template <class GridType> +void PlanarRodAssembler<GridType>:: assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol, Dune::BCRSMatrix<MatrixBlock>& matrix) { @@ -301,9 +301,9 @@ assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol, -template <class GridType, int polOrd> +template <class GridType> template <class MatrixType> -void PlanarRodAssembler<GridType, polOrd>:: +void PlanarRodAssembler<GridType>:: getLocalMatrix( EntityType &entity, const std::vector<RigidBodyMotion<2> >& localSolution, const int matSize, MatrixType& localMat) const @@ -320,7 +320,7 @@ getLocalMatrix( EntityType &entity, Dune::P1LocalFiniteElement<double,double,gridDim> localFiniteElement; // Get quadrature rule - const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(entity.type(), polOrd); + const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(entity.type(), 2); /* Loop over all integration points */ for (int ip=0; ip<quad.size(); ip++) { @@ -453,8 +453,8 @@ getLocalMatrix( EntityType &entity, } -template <class GridType, int polOrd> -void PlanarRodAssembler<GridType, polOrd>:: +template <class GridType> +void PlanarRodAssembler<GridType>:: assembleGradient(const std::vector<RigidBodyMotion<2> >& sol, Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const { @@ -483,7 +483,7 @@ assembleGradient(const std::vector<RigidBodyMotion<2> >& sol, localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)]; // Get quadrature rule - const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), polOrd); + const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), 2); for (int pt=0; pt<quad.size(); pt++) { @@ -555,8 +555,8 @@ assembleGradient(const std::vector<RigidBodyMotion<2> >& sol, } -template <class GridType, int polOrd> -double PlanarRodAssembler<GridType, polOrd>:: +template <class GridType> +double PlanarRodAssembler<GridType>:: computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const { const int maxlevel = grid_->maxLevel(); @@ -584,7 +584,7 @@ computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)]; // Get quadrature rule - const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), polOrd); + const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), 2); for (int pt=0; pt<quad.size(); pt++) { diff --git a/src/rodassembler.hh b/src/rodassembler.hh index 3f699d79ccb745ca42bd8afce5763ca387ce40c4..2cc2023c929d802f8e184cb0c74d04aabe89c6f0 100644 --- a/src/rodassembler.hh +++ b/src/rodassembler.hh @@ -79,7 +79,7 @@ public: /** \brief The FEM operator for a 2D extensible, shearable rod */ -template <class GridType, int polOrd> +template <class GridType> class PlanarRodAssembler : public GeodesicFEAssembler<typename GridType::LeafGridView, RigidBodyMotion<2> > {