From 8a1b7f872602a1aced9fcc46f5fec55de880800e Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 15 Nov 2006 10:11:13 +0000
Subject: [PATCH] further steps towards a working implementation
[[Imported from SVN: r1041]]
---
dirneucoupling.cc | 186 ++++++++++++++++++++++++++++------------------
1 file changed, 113 insertions(+), 73 deletions(-)
diff --git a/dirneucoupling.cc b/dirneucoupling.cc
index 8cb84786..399f627f 100644
--- a/dirneucoupling.cc
+++ b/dirneucoupling.cc
@@ -1,7 +1,5 @@
#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;
--
GitLab