further steps towards a working implementation

[[Imported from SVN: r1041]]

dirneucoupling.cc

 #include <config.h>
 
-//#define HAVE_IPOPT
-
 #include <dune/grid/onedgrid.hh>
 #include <dune/grid/uggrid.hh>
 
@@ -15,26 +13,22 @@
 #include <dune/common/bitfield.hh>
 #include <dune/common/configparser.hh>
 
-#include "../contact/src/contactmmgstep.hh"
+#include "../common/multigridstep.hh"
 #include "../solver/iterativesolver.hh"
 #include "../common/projectedblockgsstep.hh"
-#include "../common/geomestimator.hh"
 #include "../common/readbitfield.hh"
 #include "../common/energynorm.hh"
 #include "../common/boundarypatch.hh"
 #include "../common/prolongboundarypatch.hh"
+#include "../common/neumannassembler.hh"
 
 #include "src/quaternion.hh"
-//#include "src/rodassembler.hh"
+#include "src/rodassembler.hh"
 #include "src/configuration.hh"
 #include "src/averageinterface.hh"
 #include "src/rodsolver.hh"
 #include "src/rodwriter.hh"
 
-// Number of degrees of freedom of a correction for a rod configuration
-// 6 (x, y, z, v_1, v_2, v_3) for a spatial rod
-const int blocksize = 6;
-
 // Space dimension
 const int dim = 3;
 
@@ -44,12 +38,9 @@ using std::string;
 int main (int argc, char *argv[]) try
 {
     // Some types that I need
-    typedef BCRSMatrix<FieldMatrix<double, dim, dim> >             MatrixType;
-    typedef BlockVector<FieldVector<double, dim> >                 VectorType;
-
-    //typedef BCRSMatrix<FieldMatrix<double, blocksize, blocksize> > RodMatrixType;
-    //typedef BlockVector<FieldVector<double, blocksize> >           RodCorrectionType;
-    typedef std::vector<Configuration>                               RodSolutionType;
+    typedef BCRSMatrix<FieldMatrix<double, dim, dim> >   MatrixType;
+    typedef BlockVector<FieldVector<double, dim> >       VectorType;
+    typedef std::vector<Configuration>                   RodSolutionType;
 
     // parse data file
     ConfigParser parameterSet;
@@ -60,13 +51,15 @@ int main (int argc, char *argv[]) try
     const int maxLevel            = parameterSet.get("maxLevel", int(0));
     const int maxDirichletNeumannSteps = parameterSet.get("maxDirichletNeumannSteps", int(0));
     const int maxTrustRegionSteps = parameterSet.get("maxTrustRegionSteps", int(0));
-    const int numIt            = parameterSet.get("numIt", 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 int baseIterations   = parameterSet.get("baseIt", int(0));
+    const double mgTolerance     = parameterSet.get("tolerance", double(0));
     const double baseTolerance = parameterSet.get("baseTolerance", double(0));
+    const int initialTrustRegionRadius = parameterSet.get("initialTrustRegionRadius", int(0));
+    const double damping       = parameterSet.get("damping", double(1));
 
     // Problem settings
     std::string path = parameterSet.get("path", "xyz");
@@ -92,12 +85,7 @@ int main (int argc, char *argv[]) try
     
     AmiraMeshReader<GridType>::read(grid, path + objectName);
 
-
-    //Array<SimpleVector<BoxConstraint<dim> > > trustRegionObstacles(1);
-    //Array<BitField> hasObstacle(1);
-    Array<BitField> dirichletNodes(1);
-
-    double trustRegionRadius = 0.1;
+    std::vector<BitField> dirichletNodes(1);
 
     RodSolutionType rodX(rodGrid.size(0,1));
 
@@ -140,7 +128,7 @@ int main (int argc, char *argv[]) try
     dirichletValues[0].resize(grid.size(0, dim));
     AmiraMeshReader<int>::readFunction(dirichletValues[0], path + dirichletValuesFile);
 
-    Array<BoundaryPatch<GridType> > dirichletBoundary;
+    std::vector<BoundaryPatch<GridType> > dirichletBoundary;
     dirichletBoundary.resize(maxLevel+1);
     dirichletBoundary[0].setup(grid, 0);
     readBoundaryPatch(dirichletBoundary[0], path + dirichletNodesFile);
@@ -149,26 +137,23 @@ int main (int argc, char *argv[]) try
     dirichletNodes.resize(toplevel+1);
     for (int i=0; i<=toplevel; i++) {
         
-        dirichletNodes[i].resize( dim*grid.size(i,dim) + blocksize * rodGrid.size(i,1));
-        dirichletNodes[i].unsetAll();
+        dirichletNodes[i].resize( dim*grid.size(i,dim));
 
         for (int j=0; j<grid.size(i,dim); j++)
             for (int k=0; k<dim; k++)
                 dirichletNodes[i][j*dim+k] = dirichletBoundary[i].containsVertex(j);
         
-        for (int j=0; j<blocksize; j++)
-            dirichletNodes[i][dirichletNodes[i].size()-1-j] = true;
-        
     }
     
-    // ////////////////////////////////////////////////////////////
-    //    Create solution and rhs vectors
-    // ////////////////////////////////////////////////////////////
-    
-    VectorType totalRhs, totalCorr;
-    totalRhs.resize(grid.size(toplevel,dim) + 2*rodGrid.size(toplevel,1));
-    totalCorr.resize(grid.size(toplevel,dim) + 2*rodGrid.size(toplevel,1));
-    
+    // /////////////////////////////////////////////////////
+    //   Determine the interface boundary
+    // /////////////////////////////////////////////////////
+    std::vector<BoundaryPatch<GridType> > interfaceBoundary;
+    interfaceBoundary.resize(maxLevel+1);
+    interfaceBoundary[0].setup(grid, 0);
+    readBoundaryPatch(interfaceBoundary[0], path + interfaceNodesFile);
+    PatchProlongator<GridType>::prolong(interfaceBoundary);
+
     // //////////////////////////////////////////
     //   Assemble 3d linear elasticity problem
     // //////////////////////////////////////////
@@ -177,8 +162,11 @@ int main (int argc, char *argv[]) try
     LeafP1OperatorAssembler<GridType,double,dim> hessian3d(grid);
     hessian3d.assemble(lstiff,u,f);
 
+    // ////////////////////////////////////////////////////////////
+    //    Create solution and rhs vectors
+    // ////////////////////////////////////////////////////////////
+    
     VectorType x3d(grid.size(toplevel,dim));
-    //VectorType corr3d(grid.size(toplevel,dim));
     VectorType rhs3d(grid.size(toplevel,dim));
 
     // No external forces
@@ -191,59 +179,73 @@ int main (int argc, char *argv[]) try
             if (dirichletNodes[toplevel][i*dim+j])
                 x3d[i][j] = dirichletValues[toplevel][i][j];
 
+    // ///////////////////////////////////////////
+    //   Create a solver for the rod problem
+    // ///////////////////////////////////////////
+    RodAssembler<RodGridType> rodAssembler(rodGrid);
+    rodAssembler.setShapeAndMaterial(1, 1, 1, 2.5e5, 0.3);
+
+    RodSolver<RodGridType> rodSolver;
+    rodSolver.setup(rodGrid, 
+                    &rodAssembler,
+                    rodX,
+                    maxTrustRegionSteps,
+                    initialTrustRegionRadius,
+                    multigridIterations,
+                    mgTolerance,
+                    mu, nu1, nu2,
+                    baseIterations,
+                    baseTolerance);
+
     // ////////////////////////////////
     //   Create a multigrid solver
     // ////////////////////////////////
 
     // First create a gauss-seidel base solver
-    ProjectedBlockGSStep<MatrixType, VectorType> baseSolverStep;
+    BlockGSStep<MatrixType, VectorType> baseSolverStep;
 
     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,
+                                                       baseIterations,
+                                                       baseTolerance,
+                                                       &baseEnergyNorm,
+                                                       Solver::QUIET);
 
     // Make pre and postsmoothers
-    ProjectedBlockGSStep<MatrixType, VectorType> presmoother, postsmoother;
+    BlockGSStep<MatrixType, VectorType> presmoother, postsmoother;
 
-    ContactMMGStep<MatrixType, VectorType> contactMMGStep(totalHessian, totalCorr, totalRhs, 1);
+    MultigridStep<MatrixType, VectorType> multigridStep(*hessian3d, x3d, rhs3d, 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::QUIET;
+    multigridStep.setMGType(mu, nu1, nu2);
+    multigridStep.dirichletNodes_    = &dirichletNodes;
+    multigridStep.basesolver_        = &baseSolver;
+    multigridStep.presmoother_       = &presmoother;
+    multigridStep.postsmoother_      = &postsmoother;    
+    multigridStep.verbosity_         = Solver::REDUCED;
 
 
 
-    EnergyNorm<MatrixType, VectorType> energyNorm(contactMMGStep);
+    EnergyNorm<MatrixType, VectorType> energyNorm(multigridStep);
 
-    IterativeSolver<MatrixType, VectorType> solver;
-    solver.iterationStep = &contactMMGStep;
-    solver.numIt = numIt;
-    solver.verbosity_ = Solver::FULL;
-    solver.errorNorm_ = &energyNorm;
-    solver.tolerance_ = tolerance;
+    IterativeSolver<MatrixType, VectorType> solver(&multigridStep,
+                                                   multigridIterations,
+                                                   mgTolerance,
+                                                   &energyNorm,
+                                                   Solver::FULL);
 
     // ////////////////////////////////////
     //   Create the transfer operators
     // ////////////////////////////////////
-    for (int k=0; k<contactMMGStep.mgTransfer_.size(); k++)
-        delete(contactMMGStep.mgTransfer_[k]);
+    for (int k=0; k<multigridStep.mgTransfer_.size(); k++)
+        delete(multigridStep.mgTransfer_[k]);
     
-    contactMMGStep.mgTransfer_.resize(toplevel);
+    multigridStep.mgTransfer_.resize(toplevel);
     
-    for (int i=0; i<contactMMGStep.mgTransfer_.size(); i++){
+    for (int i=0; i<multigridStep.mgTransfer_.size(); i++){
         TruncatedMGTransfer<VectorType>* newTransferOp = new TruncatedMGTransfer<VectorType>;
         newTransferOp->setup(grid,i,i+1);
-        contactMMGStep.mgTransfer_[i] = newTransferOp;
+        multigridStep.mgTransfer_[i] = newTransferOp;
     }
 
     // /////////////////////////////////////////////////////
@@ -251,9 +253,8 @@ int main (int argc, char *argv[]) try
     // /////////////////////////////////////////////////////
 
     // Init interface value
-    Configuration lambda;
-    lambda.r = 0;
-    lambda.q = Quaternion<double>::identity();
+    Configuration referenceInterface = rodX[0];
+    Configuration lambda = referenceInterface;
 
     for (int i=0; i<maxDirichletNeumannSteps; i++) {
         
@@ -266,27 +267,63 @@ int main (int argc, char *argv[]) try
         // //////////////////////////////////////////////////
 
         rodX[0] = lambda;
-
+        rodSolver.setInitialSolution(rodX);
         rodSolver.solve();
 
+        rodX = rodSolver.getSol();
+
         // ///////////////////////////////////////////////////////////
         //   Extract Neumann values and transfer it to the 3d object
         // ///////////////////////////////////////////////////////////
+        FieldVector<double,dim> resultantForce  = rodAssembler.getResultantForce(rodX);
+        std::cout << "resultant force: " << resultantForce << std::endl;
+#if 0
+        FieldVector<double,dim> resultantTorque = rodAssembler.getResultantTorque(grid, rodX);
+#endif
+        VectorType neumannValues(grid.size(dim));
+        neumannValues = 0;
+        for (int j=0; j<neumannValues.size(); j++)
+            if (interfaceBoundary[grid.maxLevel()].containsVertex(j))
+                neumannValues[j] = resultantForce;
+
+        rhs3d = 0;
+        assembleAndAddNeumannTerm<GridType, VectorType>(interfaceBoundary[grid.maxLevel()],
+                                                        neumannValues,
+                                                        rhs3d);
 
         // ///////////////////////////////////////////////////////////
         //   Solve the Neumann problem for the 3d body
         // ///////////////////////////////////////////////////////////
+        multigridStep.setProblem(*hessian3d, x3d, rhs3d, grid.maxLevel()+1);
+        
+        solver.preprocess();
+        multigridStep.preprocess();
+        
+        solver.solve();
+        
+        x3d = multigridStep.getSol();
 
         // ///////////////////////////////////////////////////////////
         //   Extract new interface position and orientation
         // ///////////////////////////////////////////////////////////
+
         Configuration averageInterface;
-        computeAverageInterface(interface, x, averageInterface);
+//         x3d = 0;
+//         for (int i=0; i<x3d.size(); i++)
+//             x3d[i][2] = 1;
+
+        computeAverageInterface(interfaceBoundary[toplevel], x3d, averageInterface);
+
+        std::cout << "average interface: " << averageInterface << std::endl;
 
         // ///////////////////////////////////////////////////////////
         //   Compute new damped interface value
         // ///////////////////////////////////////////////////////////
 
+        for (int j=0; j<dim; j++)
+            lambda.r[j] = (1-damping) * lambda.r[j] + damping * (referenceInterface.r[j] + averageInterface.r[j]);
+        lambda.q = averageInterface.q.mult(referenceInterface.q);
+
     }
 
     // //////////////////////////////
@@ -296,6 +333,9 @@ int main (int argc, char *argv[]) try
     AmiraMeshWriter<GridType>::writeBlockVector(grid, x3d, "grid.sol");
     writeRod(rodX, "rod3d.result");
 
+    for (int i=0; i<rodX.size(); i++)
+        std::cout << rodX[i] << std::endl;
+
  } catch (Exception e) {
 
     std::cout << e << std::endl;