Skip to content
Snippets Groups Projects
Commit 1f3c7451 authored by akbib's avatar akbib Committed by akbib
Browse files

Adjust to trunk and make compile. 

This file is very old and will be removed soon when I cleaned up nbnonlincontact.cc

[[Imported from SVN: r12636]]
parent d0cb7274
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 58 additions and 83 deletions
nonlincontact.cc
+ 58
− 83
View file @ 1f3c7451
......@@ -93,13 +93,12 @@ void setTrustRegionObstacles(double trustRegionRadius,
template <class GridType, class MatrixType, class VectorType,class TrialBasis, class AnsatzBasis>
void solveTrustRegionProblem(array<OgdenAssembler<TrialBasis,AnsatzBasis>*,2>& ogdenAssembler,
NonlinearTwoBodyAssembler<GridType, VectorType>& contactAssembler,
TwoBodyAssembler<GridType, VectorType>& contactAssembler,
array<VectorType,2>& x,
array<VectorType,2>& rhs,
std::vector<VectorType>& rhs,
int level,
double trustRegionRadius,
std::vector<std::vector<BoxConstraint<field_type,dim> > >& trustRegionObstacles,
const std::vector<BitSetVector<dim> >& totalDirichletNodes,
::LoopSolver<VectorType>& solver,
MonotoneMGStep<MatrixType, VectorType>& multigridStep,
int maxTrustRegionSteps)
......@@ -117,10 +116,11 @@ void solveTrustRegionProblem(array<OgdenAssembler<TrialBasis,AnsatzBasis>*,2>& o
rhs[0] = 0;
rhs[1] = 0;
array<MatrixType,2> globalMatrix;
ogdenAssembler[0]->assembleProblem(globalMatrix[0],x[0], rhs[0]);
ogdenAssembler[1]->assembleProblem(globalMatrix[1],x[1], rhs[1]);
std::vector<const MatrixType*> globalMatrix(2);
for (int j=0; j<2; j++) {
globalMatrix[j] = new MatrixType;
ogdenAssembler[j]->assembleProblem(*const_cast<MatrixType*>(globalMatrix[j]), x[j], rhs[j]);
}
// Debug
#if 0
......@@ -134,12 +134,12 @@ void solveTrustRegionProblem(array<OgdenAssembler<TrialBasis,AnsatzBasis>*,2>& o
hessianFDCheck(x[1], *const_cast<MatrixType*>(ogdenAssembler[1].getMatrix()), ogdenAssembler[1]);
#endif
MatrixType bilinearForm;
contactAssembler.assembleJacobian(bilinearForm, &globalMatrix[0], &globalMatrix[1]);
contactAssembler.assembleJacobian(globalMatrix, bilinearForm);
//printSparseMatrix(std::cout, contactAssembler.getTransformationMatrix(), "BT", "--");
VectorType totalRhs(rhs[0].size() + rhs[1].size());
contactAssembler.assembleRightHandSide(totalRhs, rhs);
contactAssembler.assembleRightHandSide(rhs,totalRhs);
//std::cout << "totalRhs" << std::endl << totalRhs << std::endl;
......@@ -162,7 +162,7 @@ void solveTrustRegionProblem(array<OgdenAssembler<TrialBasis,AnsatzBasis>*,2>& o
setTrustRegionObstacles(trustRegionRadius,
trustRegionObstacles[level],
contactAssembler.obstacles_[level],
totalDirichletNodes[level]);
*multigridStep.ignoreNodes_);
if (dynamic_cast<MultigridStep<MatrixType,VectorType>*>(solver.iterationStep_))
dynamic_cast<MultigridStep<MatrixType,VectorType>*>(solver.iterationStep_)->setProblem(bilinearForm, totalCorr, totalRhs, level+1);
......@@ -213,7 +213,7 @@ void solveTrustRegionProblem(array<OgdenAssembler<TrialBasis,AnsatzBasis>*,2>& o
for (size_t j=0; j<2; j++) {
VectorType tmp(corr[j].size());
tmp = 0;
globalMatrix[j].umv(corr[j], tmp);
globalMatrix[j]->umv(corr[j], tmp);
modelDecrease += (rhs[j]*corr[j]) - 0.5 * (corr[j]*tmp);
}
......@@ -339,9 +339,11 @@ int main (int argc, char *argv[]) try
typedef BoundaryPatch<GridType::LevelGridView> LevelBoundaryPatch;
typedef BoundaryPatch<GridType::LeafGridView> LeafBoundaryPatch;
typedef GeometryGrid<GridType,Deformation<GridType,VectorType> > DefGridType;
typedef P1NodalBasis<GridType::LeafGridView, double> P1Basis;
typedef P1NodalBasis<DefGridType::LeafGridView, double> DefP1Basis;
typedef DeformationFunction<GridType::LeafGridView,VectorType> Deformation;
typedef GeometryGrid<GridType,Deformation> DefGridType;
typedef P1NodalBasis<GridType::LeafGridView, field_type> P1Basis;
typedef BoundaryPatch<DefGridType::LevelGridView> DefLevelBoundaryPatch;
array<GridType,2> grid;
......@@ -351,24 +353,6 @@ int main (int argc, char *argv[]) try
AmiraMeshReader<GridType>::read(grid[0], path + objectName0);
AmiraMeshReader<GridType>::read(grid[1], path + objectName1);
//from now on use the (not yet) deformed grids
array<VectorType,2> deformation;
deformation[0].resize(grid[0].size(dim));
deformation[0]=0;
deformation[1].resize(grid[1].size(dim));
deformation[1]=0;
std::vector<Deformation<GridType,VectorType>*> deform(2);
deform[0]= new Deformation<GridType,VectorType>();
deform[1]= new Deformation<GridType,VectorType>();
for (int i=0;i<2;i++){
deform[i]->setGrid(grid[i]);
deform[i]->initDeformation(deformation[i],maxLevel);
}
std::vector<DefGridType*> defGrid(2);
defGrid[0]=new DefGridType(grid[0],*deform[0]);
defGrid[1]=new DefGridType(grid[1],*deform[1]);
array<std::vector<VectorType>, 2> dirichletValues;
dirichletValues[0].resize(1);
dirichletValues[1].resize(1);
......@@ -397,13 +381,6 @@ int main (int argc, char *argv[]) try
dirichletBoundary[0][0] = coarseDirichletBoundary[0];
dirichletBoundary[1][0] = coarseDirichletBoundary[1];
// read field describing the obstacles
BitSetVector<1> obsField;
readBitField(obsField, grid[0].size(0, dim), path + obsFilename);
LevelBoundaryPatch untransfObsPatch(grid[0].levelView(0), obsField);
BoundaryPatch<DefGridType::LevelGridView> obsPatch(defGrid[0]->levelView(0), obsField);
printf("obsPatch contains %d faces\n", obsPatch.numFaces());
// //////////////////////////////////////////////////////////
// Prolong Dirichlet information to all grid levels
// //////////////////////////////////////////////////////////
......@@ -426,7 +403,7 @@ int main (int argc, char *argv[]) try
// Create solution and rhs vectors
// ////////////////////////////////////////////////////////////
array<VectorType, 2> rhs;
std::vector<VectorType> rhs(2);
array<VectorType, 2> x;
rhs[0].resize(grid[0].size(grid[0].maxLevel(), dim));
......@@ -445,32 +422,40 @@ int main (int argc, char *argv[]) try
else
AmiraMeshReader<GridType>::readFunction(x[1], interSol1);
// create deformed grids
array<Deformation*,2> deformation;
array<DefGridType*,2> defGrid;
for (int i=0;i<2;i++) {
deformation[i]=new Deformation(grid[i].leafGridView(),x[i]);
defGrid[i]= new DefGridType(&grid[i],deformation[i]);
}
// make dirichlet bitfields containing dirichlet information for both grids
std::vector<BitSetVector<dim> > totalDirichletNodes(1);
// read field describing the obstacles and the mortar boundary
BitSetVector<1> obsField;
readBitField(obsField, grid[0].size(0, dim), path + obsFilename);
LevelBoundaryPatch untransfObsPatch(grid[0].levelGridView(0), obsField);
DefLevelBoundaryPatch obsPatch(defGrid[0]->levelGridView(0), obsField);
printf("obsPatch contains %d faces\n", obsPatch.numFaces());
int offset = dirichletValues[0][0].size();
DefLevelBoundaryPatch mortarPatch(defGrid[1]->levelGridView(0),true);
if (parameterSet.hasKey("mortarPatch"))
readBoundaryPatch<DefGridType>(mortarPatch, path + parameterSet.get<string>("mortarPatch"));
totalDirichletNodes[0].resize(dirichletNodes[0][0].size() + dirichletNodes[1][0].size());
// make dirichlet bitfields containing dirichlet information for both grids
int size = grid[0].size(dim)+grid[1].size(dim);
BitSetVector<dim> totalDirichletNodes(size);
int idx = 0;
for (int j=0; j<dirichletNodes[0][0].size(); j++)
totalDirichletNodes[0][j] = dirichletNodes[0][0][j];
totalDirichletNodes[idx++] = dirichletNodes[0][0][j];
for (int j=0; j<dirichletNodes[1][0].size(); j++)
totalDirichletNodes[0][offset + j] = dirichletNodes[1][0][j];
totalDirichletNodes[idx++] = dirichletNodes[1][0][j];
// Assemble contact problem
NonlinearTwoBodyAssembler<DefGridType, VectorType> contactAssembler(*defGrid[0],
*defGrid[1],
obsPatch,
obsDistance,
CouplingPairBase::CONTACT
#ifdef EXPLICIT_CONTACT_DIRECTIONS
, &contactDirection
#endif
);
#warning The nonmortar-Lagrange matrix needs to be assembled using a quadrature rule now!
TwoBodyAssembler<DefGridType, VectorType> contactAssembler(*defGrid[0], *defGrid[1],
obsPatch, mortarPatch,
obsDistance, CouplingPairBase::CONTACT);
// //////////////////////////////////////////////////////
// Lengthen obstacle vector, to formally account for both
......@@ -478,10 +463,10 @@ int main (int argc, char *argv[]) try
// //////////////////////////////////////////////////////
std::vector<BitSetVector<1> > hasObstacle(1);
hasObstacle[0].resize(grid[0].size(0, dim) + grid[1].size(0, dim), false);
hasObstacle[0].resize(size, false);
std::vector<std::vector<BoxConstraint<field_type,dim> > > trustRegionObstacles(1);
trustRegionObstacles[0].resize(grid[0].size(0, dim) + grid[1].size(0, dim));
trustRegionObstacles[0].resize(size);
// //////////////////////////////////////
// Create a solver
......@@ -501,7 +486,7 @@ int main (int argc, char *argv[]) try
MonotoneMGStep<MatrixType, VectorType> multigridStep(1);
multigridStep.setMGType(mu, nu1, nu2);
multigridStep.ignoreNodes_ = &totalDirichletNodes[0];
multigridStep.ignoreNodes_ = &totalDirichletNodes;
multigridStep.basesolver_ = &baseSolver;
multigridStep.setSmoother(&presmoother, &postsmoother);
multigridStep.hasObstacle_ = &hasObstacle;
......@@ -580,19 +565,11 @@ int main (int argc, char *argv[]) try
//update the deformed grids!
for (int j=0;j<2;j++)
deform[j]->updateDeformation(x[j],0);
deformation[j]->setDeformation(x[j]);
contactAssembler.assembleObstacle();
contactAssembler.assembleTransferOperator();
// //////////////////////////////////////////////////////
// Lengthen obstacle vector, to formally account for both
// bodies, even though the second one doesn't actually have obstacles
// //////////////////////////////////////////////////////
for (int j=0; j<=grid[0].maxLevel(); j++)
hasObstacle[j].setAll();
/** \todo We don't really need this, because we have only one level anyways */
// //////////////////////////////////////////////////////////////////////
......@@ -603,7 +580,7 @@ int main (int argc, char *argv[]) try
ContactMGTransfer<VectorType>* newTransferOp = new ContactMGTransfer<VectorType>;
newTransferOp->setup(grid[0], grid[1], k, k+1,
contactAssembler.contactCoupling_.mortarLagrangeMatrix(k),
contactAssembler.contactCoupling_[0].mortarLagrangeMatrix(k),
contactAssembler.localCoordSystems_[k],
contactAssembler.localCoordSystems_[k+1],
hasObstacle[k], // BUG!!
......@@ -639,7 +616,7 @@ int main (int argc, char *argv[]) try
MatrixType transformedLinearMatrix;
contactAssembler.assemble(submat, transformedLinearMatrix);
dynamic_cast<EnergyNorm<MatrixType,VectorType>*>(solver.errorNorm_)->setMatrix(&transformedLinearMatrix);
energyNorm.setMatrix(&transformedLinearMatrix);
// /////////////////////////////////////////////////////
// Trust-Region Solver
......@@ -649,7 +626,6 @@ int main (int argc, char *argv[]) try
x, rhs, grid[0].maxLevel(),
trustRegionRadius,
trustRegionObstacles,
totalDirichletNodes,
solver,
multigridStep,
maxNewtonStepsI);
......@@ -728,9 +704,9 @@ int main (int argc, char *argv[]) try
p1Basis.update();
adaptor.adapt(x[i]);
deform[i]->updateDeformation(x[i],grid[i].maxLevel());
deformation[i]->setDeformation(x[i]);
deformation[i]->setGridView(grid[i].leafGridView());
defGrid[i]->update();
}
std::cout << "########################################################" << std::endl;
......@@ -814,13 +790,13 @@ int main (int argc, char *argv[]) try
int offset = dirichletValues[0][i].size();
totalDirichletNodes[i].resize(dirichletNodes[0][i].size() + dirichletNodes[1][i].size());
totalDirichletNodes.resize(dirichletNodes[0][i].size() + dirichletNodes[1][i].size());
for (int j=0; j<dirichletNodes[0][i].size(); j++)
totalDirichletNodes[i][j] = dirichletNodes[0][i][j];
totalDirichletNodes[j] = dirichletNodes[0][i][j];
for (int j=0; j<dirichletNodes[1][i].size(); j++)
totalDirichletNodes[i][offset + j] = dirichletNodes[1][i][j];
totalDirichletNodes[offset + j] = dirichletNodes[1][i][j];
}
......@@ -836,7 +812,7 @@ int main (int argc, char *argv[]) try
//update the deformed grids!
for (int j=0;j<2;j++)
deform[j]->updateDeformation(x[j],grid[j].maxLevel());
deformation[j]->setDeformation(x[j]);
contactAssembler.assembleObstacle();
contactAssembler.assembleTransferOperator();
......@@ -863,7 +839,7 @@ int main (int argc, char *argv[]) try
ContactMGTransfer<VectorType>* newTransferOp = new ContactMGTransfer<VectorType>;
newTransferOp->setup(grid[0], grid[1], k, k+1,
contactAssembler.contactCoupling_.mortarLagrangeMatrix(k),
contactAssembler.contactCoupling_[0].mortarLagrangeMatrix(k),
contactAssembler.localCoordSystems_[k],
contactAssembler.localCoordSystems_[k+1],
*contactAssembler.nonmortarBoundary_[0][k].getVertices(),
......@@ -898,7 +874,7 @@ int main (int argc, char *argv[]) try
MatrixType transformedLinearMatrix;
contactAssembler.assemble(submat, transformedLinearMatrix);
dynamic_cast<EnergyNorm<MatrixType,VectorType>*>(solver.errorNorm_)->setMatrix(&transformedLinearMatrix);
energyNorm.setMatrix(&transformedLinearMatrix);
// /////////////////////////////////////////////////////
// Trust-Region Solver
......@@ -909,7 +885,6 @@ int main (int argc, char *argv[]) try
x, rhs, grid[0].maxLevel(),
trustRegionRadius,
trustRegionObstacles,
totalDirichletNodes,
solver,
multigridStep,
maxNewtonStepsII);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment