diff --git a/dirneucoupling.cc b/dirneucoupling.cc
index 399f627fc94c1cc6b953d6b7547373ea21810745..3ac1eb753d4b5607858a142fe7e7e9773ea4387d 100644
--- a/dirneucoupling.cc
+++ b/dirneucoupling.cc
@@ -1,5 +1,7 @@
#include <config.h>
+#define HAVE_IPOPT
+
#include <dune/grid/onedgrid.hh>
#include <dune/grid/uggrid.hh>
@@ -16,6 +18,7 @@
#include "../common/multigridstep.hh"
#include "../solver/iterativesolver.hh"
#include "../common/projectedblockgsstep.hh"
+#include "../common/linearipopt.hh"
#include "../common/readbitfield.hh"
#include "../common/energynorm.hh"
#include "../common/boundarypatch.hh"
@@ -73,13 +76,13 @@ int main (int argc, char *argv[]) try
// ///////////////////////////////////////
// Create the rod grid
// ///////////////////////////////////////
- typedef OneDGrid<1,1> RodGridType;
+ typedef OneDGrid RodGridType;
RodGridType rodGrid(numRodBaseElements, 0, 5);
// ///////////////////////////////////////
// Create the grid for the 3d object
// ///////////////////////////////////////
- typedef UGGrid<dim,dim> GridType;
+ typedef UGGrid<dim> GridType;
GridType grid;
grid.setRefinementType(GridType::COPY);
@@ -202,15 +205,9 @@ int main (int argc, char *argv[]) try
// ////////////////////////////////
// First create a gauss-seidel base solver
- BlockGSStep<MatrixType, VectorType> baseSolverStep;
-
- EnergyNorm<MatrixType, VectorType> baseEnergyNorm(baseSolverStep);
-
- IterativeSolver<MatrixType, VectorType> baseSolver(&baseSolverStep,
- baseIterations,
- baseTolerance,
- &baseEnergyNorm,
- Solver::QUIET);
+ LinearIPOptSolver<VectorType> baseSolver;
+ baseSolver.verbosity_ = NumProc::QUIET;
+ baseSolver.tolerance_ = baseTolerance;
// Make pre and postsmoothers
BlockGSStep<MatrixType, VectorType> presmoother, postsmoother;
@@ -329,8 +326,8 @@ int main (int argc, char *argv[]) try
// //////////////////////////////
// Output result
// //////////////////////////////
- AmiraMeshWriter<GridType>::writeGrid(grid, "grid.result");
- AmiraMeshWriter<GridType>::writeBlockVector(grid, x3d, "grid.sol");
+ LeafAmiraMeshWriter<GridType>::writeGrid(grid, "grid.result");
+ LeafAmiraMeshWriter<GridType>::writeBlockVector(grid, x3d, "grid.sol");
writeRod(rodX, "rod3d.result");
for (int i=0; i<rodX.size(); i++)
diff --git a/rod3d.cc b/rod3d.cc
index 289f13958881939dd9ce67f2c69dcb02592a3e8d..225dea78a0a017eff24d2086794a01e8109114c4 100644
--- a/rod3d.cc
+++ b/rod3d.cc
@@ -9,8 +9,7 @@
#include "src/quaternion.hh"
#include "src/rodassembler.hh"
-#include "../common/trustregiongsstep.hh"
-#include "../contact/src/contactmmgstep.hh"
+#include "src/rodsolver.hh"
#include "../solver/iterativesolver.hh"
@@ -28,111 +27,37 @@ const int blocksize = 6;
using namespace Dune;
using std::string;
-
-void setTrustRegionObstacles(double trustRegionRadius,
- SimpleVector<BoxConstraint<blocksize> >& trustRegionObstacles)
-{
- for (int j=0; j<trustRegionObstacles.size(); j++) {
-
- for (int k=0; k<blocksize; k++) {
-
-// if (totalDirichletNodes[j*dim+k])
-// continue;
-
- trustRegionObstacles[j].val[2*k] = -trustRegionRadius;
- trustRegionObstacles[j].val[2*k+1] = trustRegionRadius;
-
- }
-
- }
-
- //std::cout << trustRegionObstacles << std::endl;
-// exit(0);
-}
-
-
int main (int argc, char *argv[]) try
{
- // Some types that I need
- typedef BCRSMatrix<FieldMatrix<double, blocksize, blocksize> > MatrixType;
- typedef BlockVector<FieldVector<double, blocksize> > CorrectionType;
- typedef std::vector<Configuration> SolutionType;
+ typedef std::vector<Configuration> SolutionType;
// parse data file
ConfigParser parameterSet;
parameterSet.parseFile("rod3d.parset");
// read solver settings
- const int minLevel = parameterSet.get("minLevel", int(0));
- const int maxLevel = parameterSet.get("maxLevel", int(0));
- const int maxNewtonSteps = parameterSet.get("maxNewtonSteps", int(0));
- const int numIt = parameterSet.get("numIt", int(0));
+ const int numLevels = parameterSet.get("numLevels", int(1));
+ const int maxTrustRegionSteps = parameterSet.get("maxNewtonSteps", int(0));
+ const int multigridIterations = parameterSet.get("numIt", int(0));
const int nu1 = parameterSet.get("nu1", int(0));
const int nu2 = parameterSet.get("nu2", int(0));
const int mu = parameterSet.get("mu", int(0));
- const int baseIt = parameterSet.get("baseIt", int(0));
- const double tolerance = parameterSet.get("tolerance", double(0));
- const double baseTolerance = parameterSet.get("baseTolerance", double(0));
-
- // Problem settings
+ const int baseIterations = parameterSet.get("baseIt", int(0));
+ const double mgTolerance = parameterSet.get("tolerance", double(0));
+ const double baseTolerance = parameterSet.get("baseTolerance", double(0));
+ const double initialTrustRegionRadius = parameterSet.get("initialTrustRegionRadius", double(1));
const int numRodBaseElements = parameterSet.get("numRodBaseElements", int(0));
// ///////////////////////////////////////
// Create the grid
// ///////////////////////////////////////
- typedef OneDGrid<1,1> GridType;
+ typedef OneDGrid GridType;
GridType grid(numRodBaseElements, 0, 1);
- grid.globalRefine(minLevel);
+ grid.globalRefine(numLevels-1);
- Array<SimpleVector<BoxConstraint<blocksize> > > trustRegionObstacles(1);
- Array<BitField> hasObstacle(1);
std::vector<BitField> dirichletNodes(1);
- // ////////////////////////////////
- // Create a multigrid solver
- // ////////////////////////////////
-
- // First create a gauss-seidel base solver
- TrustRegionGSStep<MatrixType, CorrectionType> baseSolverStep;
-
- EnergyNorm<MatrixType, CorrectionType> baseEnergyNorm(baseSolverStep);
-
- IterativeSolver<MatrixType, CorrectionType> baseSolver;
- baseSolver.iterationStep = &baseSolverStep;
- baseSolver.numIt = baseIt;
- baseSolver.verbosity_ = Solver::QUIET;
- baseSolver.errorNorm_ = &baseEnergyNorm;
- baseSolver.tolerance_ = baseTolerance;
-
- // Make pre and postsmoothers
- TrustRegionGSStep<MatrixType, CorrectionType> presmoother;
- TrustRegionGSStep<MatrixType, CorrectionType> postsmoother;
-
- ContactMMGStep<MatrixType, CorrectionType> contactMMGStep(1);
-
- contactMMGStep.setMGType(mu, nu1, nu2);
- contactMMGStep.dirichletNodes_ = &dirichletNodes;
- contactMMGStep.basesolver_ = &baseSolver;
- contactMMGStep.presmoother_ = &presmoother;
- contactMMGStep.postsmoother_ = &postsmoother;
- contactMMGStep.hasObstacle_ = &hasObstacle;
- contactMMGStep.obstacles_ = &trustRegionObstacles;
- contactMMGStep.verbosity_ = Solver::FULL;
-
-
-
- EnergyNorm<MatrixType, CorrectionType> energyNorm(contactMMGStep);
-
- IterativeSolver<MatrixType, CorrectionType> solver;
- solver.iterationStep = &contactMMGStep;
- solver.numIt = numIt;
- solver.verbosity_ = Solver::FULL;
- solver.errorNorm_ = &energyNorm;
- solver.tolerance_ = tolerance;
-
- double trustRegionRadius = 1;
-
SolutionType x(grid.size(grid.maxLevel(),1));
// //////////////////////////
@@ -169,210 +94,56 @@ int main (int argc, char *argv[]) try
//x[0].r[2] = -1;
- // /////////////////////////////////////////////////////////////////////
- // Refinement Loop
- // /////////////////////////////////////////////////////////////////////
-
- for (int toplevel=minLevel; toplevel<=maxLevel; toplevel++) {
-
- std::cout << "####################################################" << std::endl;
- std::cout << " Solving on level: " << toplevel << std::endl;
- std::cout << "####################################################" << std::endl;
-
- dirichletNodes.resize(toplevel+1);
- for (int i=0; i<=toplevel; i++) {
-
- dirichletNodes[i].resize( blocksize * grid.size(i,1), false );
-
- for (int j=0; j<blocksize; j++) {
- dirichletNodes[i][j] = true;
- dirichletNodes[i][dirichletNodes[i].size()-1-j] = true;
- }
- }
-
- // ////////////////////////////////////////////////////////////
- // Create solution and rhs vectors
- // ////////////////////////////////////////////////////////////
-
-
- MatrixType hessianMatrix;
- RodAssembler<GridType> rodAssembler(grid);
- rodAssembler.setShapeAndMaterial(0.01, 0.0001, 0.0001, 2.5e5, 0.3);
-
- std::cout << "Energy: " << rodAssembler.computeEnergy(x) << std::endl;
-
- MatrixIndexSet indices(grid.size(toplevel,1), grid.size(toplevel,1));
- rodAssembler.getNeighborsPerVertex(indices);
- indices.exportIdx(hessianMatrix);
-
- CorrectionType rhs, corr;
- rhs.resize(grid.size(toplevel,1));
- corr.resize(grid.size(toplevel,1));
-
-
- // //////////////////////////////////////////////////////////
- // Create obstacles
- // //////////////////////////////////////////////////////////
-
- hasObstacle.resize(toplevel+1);
- for (int i=0; i<hasObstacle.size(); i++) {
- hasObstacle[i].resize(grid.size(i, 1));
- hasObstacle[i].setAll();
- }
+ dirichletNodes.resize(numLevels);
+ for (int i=0; i<numLevels; i++) {
- trustRegionObstacles.resize(toplevel+1);
+ dirichletNodes[i].resize( blocksize * grid.size(i,1), false );
- for (int i=0; i<toplevel+1; i++) {
- trustRegionObstacles[i].resize(grid.size(i,1));
+ for (int j=0; j<blocksize; j++) {
+ dirichletNodes[i][j] = true;
+ dirichletNodes[i][dirichletNodes[i].size()-1-j] = true;
}
-
- trustRegionObstacles.resize(toplevel+1);
- for (int i=0; i<=toplevel; i++)
- trustRegionObstacles[i].resize(grid.size(i, 1));
-
- // ////////////////////////////////////
- // Create the transfer operators
- // ////////////////////////////////////
- for (int k=0; k<contactMMGStep.mgTransfer_.size(); k++)
- delete(contactMMGStep.mgTransfer_[k]);
-
- contactMMGStep.mgTransfer_.resize(toplevel);
-
- for (int i=0; i<contactMMGStep.mgTransfer_.size(); i++){
- TruncatedMGTransfer<CorrectionType>* newTransferOp = new TruncatedMGTransfer<CorrectionType>;
- newTransferOp->setup(grid,i,i+1);
- contactMMGStep.mgTransfer_[i] = newTransferOp;
- }
-
- // /////////////////////////////////////////////////////
- // Trust-Region Solver
- // /////////////////////////////////////////////////////
- for (int i=0; i<maxNewtonSteps; i++) {
-
- std::cout << "----------------------------------------------------" << std::endl;
- std::cout << " Trust-Region Step Number: " << i << std::endl;
- std::cout << "----------------------------------------------------" << std::endl;
-
- std::cout << "### Trust-Region Radius: " << trustRegionRadius << " ###" << std::endl;
-
- rhs = 0;
- corr = 0;
-
- rodAssembler.assembleGradient(x, rhs);
- rodAssembler.assembleMatrix(x, hessianMatrix);
-
- rhs *= -1;
-
- // Create trust-region obstacle on maxlevel
- setTrustRegionObstacles(trustRegionRadius,
- trustRegionObstacles[toplevel]);
-
- dynamic_cast<MultigridStep<MatrixType,CorrectionType>*>(solver.iterationStep)->setProblem(hessianMatrix, corr, rhs, toplevel+1);
-
- solver.preprocess();
-
- contactMMGStep.preprocess();
-
-
- // /////////////////////////////
- // Solve !
- // /////////////////////////////
- solver.solve();
-
- corr = contactMMGStep.getSol();
-
- //std::cout << "Correction: " << std::endl << corr << std::endl;
-
- printf("infinity norm of the correction: %g\n", corr.infinity_norm());
- if (corr.infinity_norm() < 1e-5) {
- std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
- break;
- }
-
- // ////////////////////////////////////////////////////
- // Check whether trust-region step can be accepted
- // ////////////////////////////////////////////////////
-
- SolutionType newIterate = x;
- for (int j=0; j<newIterate.size(); j++) {
-
- // Add translational correction
- for (int k=0; k<3; k++)
- newIterate[j].r[k] += corr[j][k];
-
- // Add rotational correction
- Quaternion<double> qCorr = Quaternion<double>::exp(corr[j][3], corr[j][4], corr[j][5]);
- newIterate[j].q = newIterate[j].q.mult(qCorr);
-
- }
-
-#if 0
- std::cout << "newIterate: \n";
- for (int j=0; j<newIterate.size(); j++)
- printf("%d: (%g %g %g) (%g %g %g %g)\n", j,
- newIterate[j].r[0],newIterate[j].r[1],newIterate[j].r[2],
- newIterate[j].q[0],newIterate[j].q[1],newIterate[j].q[2],newIterate[j].q[3]);
-#endif
-
- /** \todo Don't always recompute oldEnergy */
- double oldEnergy = rodAssembler.computeEnergy(x);
- double energy = rodAssembler.computeEnergy(newIterate);
-
- // compute the model decrease
- // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
- // Note that rhs = -g
- CorrectionType tmp(corr.size());
- tmp = 0;
- hessianMatrix.mmv(corr, tmp);
- double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
-
- std::cout << "Model decrease: " << modelDecrease
- << ", functional decrease: " << oldEnergy - energy << std::endl;
-
- assert(modelDecrease >= 0);
-
- if (energy >= oldEnergy) {
- printf("Richtung ist keine Abstiegsrichtung!\n");
-// std::cout << "corr[0]\n" << corr[0] << std::endl;
- //exit(0);
- }
-
- // //////////////////////////////////////////////
- // Check for acceptance of the step
- // //////////////////////////////////////////////
- if ( (oldEnergy-energy) / modelDecrease > 0.9) {
- // very successful iteration
-
- x = newIterate;
- trustRegionRadius *= 2;
-
- } else if ( (oldEnergy-energy) / modelDecrease > 0.01) {
- // successful iteration
- x = newIterate;
-
- } else {
- // unsuccessful iteration
- trustRegionRadius /= 2;
- std::cout << "Unsuccessful iteration!" << std::endl;
- }
-
- // Write current energy
- std::cout << "--- Current energy: " << energy << " ---" << std::endl;
- }
-
- // //////////////////////////////
- // Output result
- // //////////////////////////////
- writeRod(x, "rod3d.result");
- BlockVector<FieldVector<double, 6> > strain(x.size()-1);
- rodAssembler.getStrain(x,strain);
- //std::cout << strain << std::endl;
- //exit(0);
-
- writeRod(x, strain, "rod3d.strain");
-
}
+
+ // ///////////////////////////////////////////
+ // Create a solver for the rod problem
+ // ///////////////////////////////////////////
+ RodAssembler<GridType> rodAssembler(grid);
+ rodAssembler.setShapeAndMaterial(0.01, 0.0001, 0.0001, 2.5e5, 0.3);
+
+ RodSolver<GridType> rodSolver;
+ rodSolver.setup(grid,
+ &rodAssembler,
+ x,
+ maxTrustRegionSteps,
+ initialTrustRegionRadius,
+ multigridIterations,
+ mgTolerance,
+ mu, nu1, nu2,
+ baseIterations,
+ baseTolerance);
+
+ // /////////////////////////////////////////////////////
+ // Solve!
+ // /////////////////////////////////////////////////////
+
+ std::cout << "Energy: " << rodAssembler.computeEnergy(x) << std::endl;
+
+ rodSolver.setInitialSolution(x);
+ rodSolver.solve();
+ x = rodSolver.getSol();
+
+ // //////////////////////////////
+ // Output result
+ // //////////////////////////////
+ writeRod(x, "rod3d.result");
+ BlockVector<FieldVector<double, 6> > strain(x.size()-1);
+ rodAssembler.getStrain(x,strain);
+ //std::cout << strain << std::endl;
+ //exit(0);
+
+ writeRod(x, strain, "rod3d.strain");
} catch (Exception e) {
diff --git a/staticrod.cc b/staticrod.cc
index d5f472fae351c1ad735aceeca919581fea8a3ef7..46190399f1cd14d1402ec5949d60f90af787a1ab 100644
--- a/staticrod.cc
+++ b/staticrod.cc
@@ -28,8 +28,8 @@ using namespace Dune;
using std::string;
void setTrustRegionObstacles(double trustRegionRadius,
- SimpleVector<BoxConstraint<blocksize> >& trustRegionObstacles,
- const SimpleVector<BoxConstraint<blocksize> >& trueObstacles,
+ std::vector<BoxConstraint<blocksize> >& trustRegionObstacles,
+ const std::vector<BoxConstraint<blocksize> >& trueObstacles,
const BitField& dirichletNodes)
{
//std::cout << "True obstacles\n" << trueObstacles << std::endl;
@@ -86,7 +86,7 @@ int main (int argc, char *argv[]) try
// ///////////////////////////////////////
// Create the two grids
// ///////////////////////////////////////
- typedef OneDGrid<1,1> RodGridType;
+ typedef OneDGrid RodGridType;
RodGridType rod(numRodBaseElements, 0, 1);
// refine uniformly until maxLevel
@@ -97,7 +97,7 @@ int main (int argc, char *argv[]) try
int numRodElements = rod.size(maxlevel, 0);
- Array<BitField> dirichletNodes;
+ std::vector<BitField> dirichletNodes;
dirichletNodes.resize(maxLevel+1);
for (int i=0; i<=maxlevel; i++) {
@@ -156,8 +156,8 @@ int main (int argc, char *argv[]) try
hasObstacle[i].setAll();
}
- Array<SimpleVector<BoxConstraint<3> > > trueObstacles(maxlevel+1);
- Array<SimpleVector<BoxConstraint<3> > > trustRegionObstacles(maxlevel+1);
+ Array<std::vector<BoxConstraint<3> > > trueObstacles(maxlevel+1);
+ Array<std::vector<BoxConstraint<3> > > trustRegionObstacles(maxlevel+1);
for (int i=0; i<maxlevel+1; i++) {
trueObstacles[i].resize(rod.size(i,1));
@@ -179,12 +179,11 @@ int main (int argc, char *argv[]) try
EnergyNorm<MatrixType, VectorType> baseEnergyNorm(baseSolverStep);
- IterativeSolver<MatrixType, VectorType> baseSolver;
- baseSolver.iterationStep = &baseSolverStep;
- baseSolver.numIt = baseIt;
- baseSolver.verbosity_ = Solver::QUIET;
- baseSolver.errorNorm_ = &baseEnergyNorm;
- baseSolver.tolerance_ = baseTolerance;
+ IterativeSolver<MatrixType, VectorType> baseSolver(&baseSolverStep,
+ baseIt,
+ baseTolerance,
+ &baseEnergyNorm,
+ Solver::QUIET);
// Make pre and postsmoothers
ProjectedBlockGSStep<MatrixType, VectorType> presmoother;
@@ -210,12 +209,11 @@ int main (int argc, char *argv[]) try
EnergyNorm<MatrixType, VectorType> energyNorm(contactMMGStep);
- IterativeSolver<MatrixType, VectorType> solver;
- solver.iterationStep = &contactMMGStep;
- solver.numIt = numIt;
- solver.verbosity_ = Solver::QUIET;
- solver.errorNorm_ = &energyNorm;
- solver.tolerance_ = tolerance;
+ IterativeSolver<MatrixType, VectorType> solver(&contactMMGStep,
+ numIt,
+ tolerance,
+ &energyNorm,
+ Solver::QUIET);
// ///////////////////////////////////////////////////
// Do a homotopy of the material parameters
@@ -267,7 +265,7 @@ int main (int argc, char *argv[]) try
//std::cout << "Trust Region obstacles:" << std::endl;
//std::cout << (*contactMMGStep.obstacles_)[maxlevel] << std::endl;
- solver.iterationStep->setProblem(hessianMatrix, corr, rhs);
+ solver.iterationStep_->setProblem(hessianMatrix, corr, rhs);
solver.preprocess();
contactMMGStep.preprocess();
diff --git a/staticrod2.cc b/staticrod2.cc
index 48dbbce8c8a9c0fd3dea9ca94529d6b44bbe2421..0143b22c8b75e9cc645eea73f6cce06c0e709ebd 100644
--- a/staticrod2.cc
+++ b/staticrod2.cc
@@ -1,21 +1,23 @@
#include <config.h>
//#define DUNE_EXPRESSIONTEMPLATES
+#include <dune/common/bitfield.hh>
+
#include <dune/grid/onedgrid.hh>
#include <dune/istl/io.hh>
+#include <dune/disc/operators/p1operator.hh>
+
#include "../common/boundarypatch.hh"
-#include <dune/common/bitfield.hh>
-#include "src/rodassembler.hh"
+#include "src/planarrodassembler.hh"
#include "../common/projectedblockgsstep.hh"
#include "../contact/src/contactmmgstep.hh"
-#include <dune/solver/iterativesolver.hh>
+#include "../solver/iterativesolver.hh"
#include "../common/geomestimator.hh"
-#include "../common/refinegrid.hh"
#include "../common/energynorm.hh"
#include "src/rodwriter.hh"
@@ -29,7 +31,7 @@ using namespace Dune;
using std::string;
void setTrustRegionObstacles(double trustRegionRadius,
- SimpleVector<BoxConstraint<blocksize> >& trustRegionObstacles,
+ std::vector<BoxConstraint<blocksize> >& trustRegionObstacles,
const SimpleVector<BoxConstraint<blocksize> >& trueObstacles,
const BitField& dirichletNodes)
{
@@ -95,12 +97,12 @@ int main (int argc, char *argv[]) try
// ///////////////////////////////////////
// Create the two grids
// ///////////////////////////////////////
- typedef OneDGrid<1,1> GridType;
+ typedef OneDGrid GridType;
GridType grid(numRodBaseElements, 0, 1);
- Array<SimpleVector<BoxConstraint<3> > > trustRegionObstacles(1);
+ Array<std::vector<BoxConstraint<3> > > trustRegionObstacles(1);
Array<BitField> hasObstacle(1);
- Array<BitField> dirichletNodes(1);
+ std::vector<BitField> dirichletNodes(1);
// ////////////////////////////////
// Create a multigrid solver
@@ -111,18 +113,17 @@ int main (int argc, char *argv[]) try
EnergyNorm<MatrixType, VectorType> baseEnergyNorm(baseSolverStep);
- IterativeSolver<MatrixType, VectorType> baseSolver;
- baseSolver.iterationStep = &baseSolverStep;
- baseSolver.numIt = baseIt;
- baseSolver.verbosity_ = Solver::QUIET;
- baseSolver.errorNorm_ = &baseEnergyNorm;
- baseSolver.tolerance_ = baseTolerance;
+ IterativeSolver<MatrixType, VectorType> baseSolver(&baseSolverStep,
+ baseIt,
+ baseTolerance,
+ &baseEnergyNorm,
+ Solver::QUIET);
// Make pre and postsmoothers
ProjectedBlockGSStep<MatrixType, VectorType> presmoother;
ProjectedBlockGSStep<MatrixType, VectorType> postsmoother;
- ContactMMGStep<MatrixType, VectorType, GridType > contactMMGStep(1);
+ ContactMMGStep<MatrixType, VectorType> contactMMGStep(1);
contactMMGStep.setMGType(mu, nu1, nu2);
contactMMGStep.dirichletNodes_ = &dirichletNodes;
@@ -137,12 +138,11 @@ int main (int argc, char *argv[]) try
EnergyNorm<MatrixType, VectorType> energyNorm(contactMMGStep);
- IterativeSolver<MatrixType, VectorType> solver;
- solver.iterationStep = &contactMMGStep;
- solver.numIt = numIt;
- solver.verbosity_ = Solver::FULL;
- solver.errorNorm_ = &energyNorm;
- solver.tolerance_ = tolerance;
+ IterativeSolver<MatrixType, VectorType> solver(&contactMMGStep,
+ numIt,
+ tolerance,
+ &energyNorm,
+ Solver::FULL);
double trustRegionRadius = 0.1;
@@ -184,7 +184,7 @@ int main (int argc, char *argv[]) try
MatrixType hessianMatrix;
- RodAssembler<GridType,4> rodAssembler(grid);
+ PlanarRodAssembler<GridType,4> rodAssembler(grid);
rodAssembler.setParameters(1, 350000, 350000);
@@ -263,7 +263,7 @@ int main (int argc, char *argv[]) try
trueObstacles[toplevel],
dirichletNodes[toplevel]);
- dynamic_cast<MultigridStep<MatrixType,VectorType,GridType>*>(solver.iterationStep)->setProblem(hessianMatrix, corr, rhs, toplevel+1);
+ dynamic_cast<MultigridStep<MatrixType,VectorType>*>(solver.iterationStep_)->setProblem(hessianMatrix, corr, rhs, toplevel+1);
solver.preprocess();
@@ -334,8 +334,20 @@ int main (int argc, char *argv[]) try
GeometricEstimator<GridType> estimator;
estimator.estimate(grid, (toplevel<=minLevel) ? refineAll : refineCondition);
- refineGridAndTransferFunction(grid, x);
-
+
+ P1FunctionManager<GridType,double> functionManager(grid);
+ LeafP1Function<GridType,double,blocksize> sol(grid);
+ *sol = x;
+
+ grid.preAdapt();
+ sol.preAdapt();
+ grid.adapt();
+
+ sol.postAdapt(functionManager);
+ grid.postAdapt();
+
+ x = *sol;
+
//writeRod(x, "solutions/rod_1.result");
}