Skip to content
Snippets Groups Projects
Commit 4f5ee53e authored by Oliver Sander's avatar Oliver Sander Committed by sander@FU-BERLIN.DE
Browse files

use the function space basis instead of the grid index set to get global indices 

[[Imported from SVN: r7867]]
parent 8ff64810
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gfe/cosseratenergystiffness.hh
+ 36
− 3
View file @ 4f5ee53e
......@@ -82,6 +82,8 @@ public: // for testing
// This second derivative is almost given by the method getFirstDerivativesOfDirectors.
// However, since the directors of a given unit quaternion are the _columns_ of the
// corresponding orthogonal matrix, we need to invert the i and j indices
//
// So, if I am not mistaken, DR[i][j][k] contains \partial M_ij / \partial k
Tensor3<double,3 , 3, 4> dd_dq;
value.q.getFirstDerivativesOfDirectors(dd_dq);
......@@ -136,7 +138,8 @@ public:
RT curvatureEnergy(const Tensor3<double,3,3,3>& DR) const
{
return mu_ * std::pow(L_c_ * curl(DR).frobenius_norm(),q_);
//return mu_ * std::pow(L_c_ * curl(DR).frobenius_norm(),q_);
return mu_ * std::pow(L_c_ * DR.frobenius_norm(),q_);
}
RT bendingEnergy(const Dune::FieldMatrix<double,dim,dim>& R, const Tensor3<double,3,3,3>& DR) const
......@@ -180,7 +183,7 @@ public:
/** \brief Curvature exponent */
double q_;
const BoundaryPatch<GridView>* neumannBoundary_;
};
template <class GridView, class LocalFiniteElement, int dim>
......@@ -195,7 +198,7 @@ energy(const Entity& element,
LocalGeodesicFEFunction<gridDim, double, LocalFiniteElement, TargetSpace> localGeodesicFEFunction(localFiniteElement,
localSolution);
int quadOrder = 1;//gridDim;
int quadOrder = 2*gridDim;
const Dune::QuadratureRule<double, gridDim>& quad
= Dune::QuadratureRules<double, gridDim>::rule(element.type(), quadOrder);
......@@ -276,7 +279,37 @@ energy(const Entity& element,
DUNE_THROW(Dune::NotImplemented, "CosseratEnergyStiffness for 1d grids");
}
//////////////////////////////////////////////////////////////////////////////
// Assemble boundary contributions
//////////////////////////////////////////////////////////////////////////////
for (typename Entity::LeafIntersectionIterator it = element.ileafbegin(); it != element.ileafend(); ++it) {
if (not neumannBoundary_ or not neumannBoundary_->contains(*it))
continue;
const Dune::QuadratureRule<double, gridDim-1>& quad
= Dune::QuadratureRules<double, gridDim-1>::rule(it->type(), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
// Local position of the quadrature point
const Dune::FieldVector<double,gridDim>& quadPos = it->geometryInInside().global(quad[pt].position());
const double integrationElement = it->geometry().integrationElement(quad[pt].position());
// The value of the local function
RigidBodyMotion<dim> value = localGeodesicFEFunction.evaluate(quadPos);
// Constant Neumann force in z-direction
energy += thickness_ * 40 * value.r[2] * quad[pt].weight() * integrationElement;
}
}
return energy;
}
......
dune/gfe/riemanniantrsolver.cc
+ 27
− 6
View file @ 4f5ee53e
......@@ -62,7 +62,9 @@ setup(const GridType& grid,
// Create a multigrid solver
// ////////////////////////////////
#if 0
// First create a Gauss-seidel base solver
TrustRegionGSStep<MatrixType, CorrectionType>* baseSolverStep = new TrustRegionGSStep<MatrixType, CorrectionType>;
EnergyNorm<MatrixType, CorrectionType>* baseEnergyNorm = new EnergyNorm<MatrixType, CorrectionType>(*baseSolverStep);
......@@ -72,6 +74,11 @@ setup(const GridType& grid,
baseTolerance,
baseEnergyNorm,
Solver::QUIET);
#else
QuadraticIPOptSolver<MatrixType,CorrectionType>* baseSolver = new QuadraticIPOptSolver<MatrixType,CorrectionType>();
baseSolver->tolerance_ = baseTolerance;
baseSolver->verbosity_ = Solver::QUIET;
#endif
// Make pre and postsmoothers
TrustRegionGSStep<MatrixType, CorrectionType>* presmoother = new TrustRegionGSStep<MatrixType, CorrectionType>;
......@@ -105,12 +112,21 @@ setup(const GridType& grid,
h1SemiNorm_ = new H1SemiNorm<CorrectionType>(*A);
#if 0
innerSolver_ = std::shared_ptr<LoopSolver<CorrectionType> >(new ::LoopSolver<CorrectionType>(mmgStep,
innerIterations_,
innerTolerance_,
h1SemiNorm_,
Solver::FULL));
#else
QuadraticIPOptSolver<MatrixType,CorrectionType>* solver = new QuadraticIPOptSolver<MatrixType,CorrectionType>();
solver->tolerance_ = innerTolerance_;
solver->verbosity_ = Solver::REDUCED;
innerSolver_ = std::shared_ptr<QuadraticIPOptSolver<MatrixType,CorrectionType> >(solver);
#endif
// Write all intermediate solutions, if requested
if (instrumented_
&& dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
......@@ -190,7 +206,7 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
std::vector<std::vector<BoxConstraint<field_type,blocksize> > > trustRegionObstacles((mgStep)
? mgStep->numLevels_
: 0);
: 1);
// /////////////////////////////////////////////////////
// Set up the log file, if requested
......@@ -242,8 +258,6 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
// The right hand side is the _negative_ gradient
rhs *= -1;
/* std::cout << "rhs:\n" << rhs << std::endl;
std::cout << "matrix[0][0]:\n" << (*hessianMatrix_)[0][0] << std::endl;*/
// //////////////////////////////////////////////////////////////////////
// Modify matrix and right-hand side to account for Dirichlet values
......@@ -277,10 +291,17 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
}
mgStep->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1);
//std::cout << "rhs:\n" << rhs << std::endl;
//std::cout << "matrix[0][0]:\n" << (*hessianMatrix_)[0][0] << std::endl;
//printmatrix(std::cout, (*hessianMatrix_), "hessian", "--");
//writeMatrixToMatlab(*hessianMatrix_, "cosserat_hessian");
//mgStep->setProblem(*hessianMatrix_, corr, rhs, grid_->maxLevel()+1);
std::dynamic_pointer_cast<QuadraticIPOptSolver<MatrixType,CorrectionType> >(innerSolver_)->setProblem(*hessianMatrix_, corr, rhs);
trustRegionObstacles.back() = trustRegion.obstacles();
mgStep->obstacles_ = &trustRegionObstacles;
//mgStep->obstacles_ = &trustRegionObstacles;
std::dynamic_pointer_cast<QuadraticIPOptSolver<MatrixType,CorrectionType> >(innerSolver_)->obstacles_ = &trustRegionObstacles.back();
innerSolver_->preprocess();
......@@ -420,7 +441,7 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
}
if (energy >= oldEnergy &&
(std::abs(oldEnergy-energy)/energy < 1e-9 || modelDecrease/energy < 1e-9)) {
(std::abs((oldEnergy-energy)/energy) < 1e-9 || std::abs(modelDecrease/energy) < 1e-9)) {
if (this->verbosity_ == NumProc::FULL)
std::cout << "Suspecting rounding problems" << std::endl;
......
dune/gfe/rodassembler.cc
+ 21
− 31
View file @ 4f5ee53e
......@@ -18,9 +18,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
{
using namespace Dune;
const typename GridView::Traits::IndexSet& indexSet = this->basis_.getGridView().indexSet();
if (sol.size()!=indexSet.size(gridDim))
if (sol.size()!=this->basis_.size())
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
grad.resize(sol.size());
......@@ -39,7 +37,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
std::vector<RigidBodyMotion<3> > localSolution(nDofs);
for (int i=0; i<nDofs; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
localSolution[i] = sol[this->basis_.index(*it,i)];
// Assemble local gradient
std::vector<FieldVector<double,blocksize> > localGradient(nDofs);
......@@ -51,7 +49,7 @@ assembleGradient(const std::vector<RigidBodyMotion<3> >& sol,
// Add to global gradient
for (int i=0; i<nDofs; i++)
grad[indexSet.subIndex(*it,i,gridDim)] += localGradient[i];
grad[this->basis_.index(*it,i)] += localGradient[i];
}
......@@ -65,17 +63,17 @@ getStrain(const std::vector<RigidBodyMotion<3> >& sol,
{
using namespace Dune;
const typename GridView::Traits::IndexSet& indexSet = this->gridView_.indexSet();
const typename GridView::Traits::IndexSet& indexSet = this->basis_.getGridView().indexSet();
if (sol.size()!=indexSet.size(gridDim))
if (sol.size()!=this->basis_.size())
DUNE_THROW(Exception, "Solution vector doesn't match the grid!");
// Strain defined on each element
strain.resize(indexSet.size(0));
strain = 0;
ElementIterator it = this->gridView_.template begin<0>();
ElementIterator endIt = this->gridView_.template end<0>();
ElementIterator it = this->basis_.getGridView().template begin<0>();
ElementIterator endIt = this->basis_.getGridView().template end<0>();
// Loop over all elements
for (; it!=endIt; ++it) {
......@@ -227,15 +225,13 @@ void RodAssembler<GridView,2>::
assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
Dune::BCRSMatrix<MatrixBlock>& matrix)
{
const typename GridView::IndexSet& indexSet = this->gridView_.indexSet();
Dune::MatrixIndexSet neighborsPerVertex;
this->getNeighborsPerVertex(neighborsPerVertex);
matrix = 0;
ElementIterator it = this->gridView_.template begin<0>();
ElementIterator endit = this->gridView_.template end<0> ();
ElementIterator it = this->basis_.getGridView().template begin<0>();
ElementIterator endit = this->basis_.getGridView().template end<0> ();
Dune::Matrix<MatrixBlock> mat;
......@@ -247,7 +243,7 @@ assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
std::vector<RigidBodyMotion<2> > localSolution(numOfBaseFct);
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
localSolution[i] = sol[this->basis_.index(*it,i)];
// setup matrix
getLocalMatrix( *it, localSolution, mat);
......@@ -255,11 +251,11 @@ assembleMatrix(const std::vector<RigidBodyMotion<2> >& sol,
// Add element matrix to global stiffness matrix
for(int i=0; i<numOfBaseFct; i++) {
int row = indexSet.subIndex(*it,i,gridDim);
int row = this->basis_.index(*it,i);
for (int j=0; j<numOfBaseFct; j++ ) {
int col = indexSet.subIndex(*it,j,gridDim);
int col = this->basis_.index(*it,j);
matrix[row][col] += mat[i][j];
}
......@@ -280,8 +276,6 @@ getLocalMatrix( EntityType &entity,
const std::vector<RigidBodyMotion<2> >& localSolution,
Dune::Matrix<MatrixBlock>& localMat) const
{
const typename GridView::IndexSet& indexSet = this->gridView_.indexSet();
/* ndof is the number of vectors of the element */
int ndof = localSolution.size();
......@@ -429,16 +423,14 @@ void RodAssembler<GridView,2>::
assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
Dune::BlockVector<Dune::FieldVector<double, blocksize> >& grad) const
{
const typename GridView::IndexSet& indexSet = this->gridView_.indexSet();
if (sol.size()!=this->gridView_.size(gridDim))
if (sol.size()!=this->basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
grad.resize(sol.size());
grad = 0;
ElementIterator it = this->gridView_.template begin<0>();
ElementIterator endIt = this->gridView_.template end<0>();
ElementIterator it = this->basis_.getGridView().template begin<0>();
ElementIterator endIt = this->basis_.getGridView().template end<0>();
// Loop over all elements
for (; it!=endIt; ++it) {
......@@ -450,7 +442,7 @@ assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
RigidBodyMotion<2> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
localSolution[i] = sol[this->basis_.index(*it,i)];
// Get quadrature rule
const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), 2);
......@@ -500,7 +492,7 @@ assembleGradient(const std::vector<RigidBodyMotion<2> >& sol,
for (int dof=0; dof<numOfBaseFct; dof++) {
int globalDof = indexSet.subIndex(*it,dof,gridDim);
int globalDof = this->basis_.index(*it,dof);
//printf("globalDof: %d partA1: %g partA3: %g\n", globalDof, partA1, partA3);
......@@ -531,13 +523,11 @@ computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const
{
double energy = 0;
const typename GridView::IndexSet& indexSet = this->gridView_.indexSet();
if (sol.size()!=this->gridView_.size(gridDim))
if (sol.size()!=this->basis_.size())
DUNE_THROW(Dune::Exception, "Solution vector doesn't match the grid!");
ElementIterator it = this->gridView_.template begin<0>();
ElementIterator endIt = this->gridView_.template end<0>();
ElementIterator it = this->basis_.getGridView().template begin<0>();
ElementIterator endIt = this->basis_.getGridView().template end<0>();
// Loop over all elements
for (; it!=endIt; ++it) {
......@@ -550,7 +540,7 @@ computeEnergy(const std::vector<RigidBodyMotion<2> >& sol) const
RigidBodyMotion<2> localSolution[numOfBaseFct];
for (int i=0; i<numOfBaseFct; i++)
localSolution[i] = sol[indexSet.subIndex(*it,i,gridDim)];
localSolution[i] = sol[this->basis_.index(*it,i)];
// Get quadrature rule
const Dune::QuadratureRule<double, gridDim>& quad = Dune::QuadratureRules<double, gridDim>::rule(it->type(), 2);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment