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

Add new method which uses the new "Problem" classes to assemble the defect...

Add new method which uses the new  "Problem" classes to assemble the defect problem. This way one can also use the estimator for problems where the energy functional is composed of multiple terms, as for example in the dynamic case.

This method is not perfect yet, as one might want to use it with ConformingBasis as well.

[[Imported from SVN: r12675]]
parent 506dbd60
No related branches found
No related tags found
No related merge requests found
......@@ -8,6 +8,9 @@
#include <dune/fufem/functionspacebases/p2hierarchicalbasis.hh>
#include <dune/fufem/boundarypatchprolongator.hh>
#include <dune/fufem/estimators/refinementindicator.hh>
#include <dune/fufem/assemblers/functionalassembler.hh>
#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
#include <dune/contact/assemblers/leafp2mortarcoupling.hh>
#include <dune/contact/assemblers/p2onebodyassembler.hh>
......@@ -361,6 +364,397 @@ void HierarchicContactEstimator<GridType,field_type>::estimate(const std::vector
delete(coupling[i]);
}
template <class GridType, class field_type>
template <class ProblemType>
void HierarchicContactEstimator<GridType,field_type>::estimate(const std::vector<VectorType>& x,
const std::vector<P1BasisGridFunction*>& volumeTerm,
const std::vector<P1BasisGridFunction*>& neumannTerm,
const std::vector<LeafBoundaryPatch>& dirichletBoundary,
const std::vector<LeafBoundaryPatch>* neumannBoundary,
std::vector<RefinementIndicator<GridType>* >& refinementIndicator,
ProblemType& contactProblem)
{
assert(x.size() == volumeTerm.size());
assert(x.size() == neumannTerm.size());
// ////////////////////////////////////////////////////////////////////
// embed x into a second order fe space with the nodal basis
// ////////////////////////////////////////////////////////////////////
for (int i=0; i<numGrids(); i++)
assert(x[i].size() == grids_[i]->leafIndexSet().size(dim));
std::vector<P2Basis*> p2NodalBasis(numGrids());
std::vector<VectorType> xP2(p2NodalBasis.size());
for (int i=0; i<numGrids(); i++) {
P1Basis p1NodalBasis(grids_[i]->leafView());
p2NodalBasis[i] = new P2Basis(grids_[i]->leafView());
xP2[i].resize(p2NodalBasis[i]->size());
P1BasisGridFunction xFunction(p1NodalBasis, x[i]);
Functions::interpolate(*p2NodalBasis[i], xP2[i], xFunction);
}
// ////////////////////////////////////////////////////////////////////
// Assemble the separate preconditioned defect problems
// ////////////////////////////////////////////////////////////////////
// ///////////////////////////////////////////////////////////////////
// Assemble external forces
// /////////////////////////////////////////////////////////////////
std::vector<VectorType> residual(numGrids());
std::vector<DiagonalMatrixType> matrix(numGrids());
Dune::Timer watch;
// assemble external forces
std::vector<VectorType> externalForces(numGrids());
for (int i=0; i<numGrids(); i++) {
externalForces[i].resize(xP2[i].size());
externalForces[i] = 0;
if (volumeTerm[i]) {
// local assembler for the volume term
L2FunctionalAssembler<GridType, typename P2Basis::LocalFiniteElement, Dune::FieldVector<field_type,dim> > localL2Assembler(*volumeTerm[i]);
FunctionalAssembler<P2Basis> funcAssembler(*p2NodalBasis[i]);
funcAssembler.assemble(localL2Assembler,externalForces[i]);
}
if (neumannTerm[i]) {
BoundaryFunctionalAssembler<P2Basis> boundaryAssembler(*p2NodalBasis[i],(*neumannBoundary)[i]);
NeumannBoundaryAssembler<GridType, Dune::FieldVector<field_type,dim> > localNeumannAssembler(*neumannTerm[i]);
boundaryAssembler.assemble(localNeumannAssembler, externalForces[i], false);
}
}
contactProblem.setExternalForces(externalForces);
contactProblem.assembleDefectQP(xP2,residual,matrix);
std::cout << watch.elapsed() << " to assemble the preconditioned defect problem." << std::endl;
// ///////////////////////////////////////////////////////////
// Construct the set of Dirichlet dofs for the p2 problems
// ///////////////////////////////////////////////////////////
std::vector<Dune::BitSetVector<dim> > p2DirichletNodes(numGrids());
for (int i=0; i<numGrids(); i++)
constructBoundaryDofs(dirichletBoundary[i], *p2NodalBasis[i], p2DirichletNodes[i]);
// //////////////////////////////////////////////////////////////////////
// Set up solution vectors for the preconditioned P2 defect problem
// //////////////////////////////////////////////////////////////////////
predXNodal_.resize(numGrids());
for (int i=0; i<numGrids(); i++) {
predXNodal_[i].resize(xP2[i].size());
predXNodal_[i] = 0;
}
// ///////////////////////////////////////////////////////
// Construct the P2 mortar coupling matrices
// ///////////////////////////////////////////////////////
// Not nice: Get the deformed grids from the contact assembler inside the problem
typedef typename ProblemType::DeformedGridType DeformedGridType;
typename ProblemType::ContactAssembler* contactAssembler = contactProblem.contactAssembler();
std::vector<LeafP2MortarCoupling<field_type,DeformedGridType>*> p2couplings(numCouplings());
for (int i=0; i<numCouplings(); i++) {
CouplingPair<DeformedGridType,DeformedGridType,field_type>& coupling = contactAssembler->coupling_[i];
if (coupling.type_ != CouplingPairBase::RIGID) {
int grid0Idx = coupling.gridIdx_[0];
int grid1Idx = coupling.gridIdx_[1];
BoundaryPatch<typename DeformedGridType::LeafGridView> leafNonmortar;
BoundaryPatchProlongator<DeformedGridType>::prolong(*coupling.patch0(),leafNonmortar);
BoundaryPatch<typename DeformedGridType::LeafGridView> leafMortar;
BoundaryPatchProlongator<DeformedGridType>::prolong(*coupling.patch1(),leafMortar);
// Create zero deformation as we have already deformed the underlying grid, this is very stupid..
typedef P2NodalBasis<typename DeformedGridType::LeafGridView> DefP2Basis;
DefP2Basis defP2Basis0(contactAssembler->grids_[grid0Idx]->leafGridView());
VectorType dummy(defP2Basis0.size());
dummy = 0;
BasisGridFunction<DefP2Basis,VectorType> def0(Functions::makeFunction(defP2Basis0,dummy));
DefP2Basis defP2Basis1(contactAssembler->grids_[grid1Idx]->leafGridView());
dummy.resize(defP2Basis1.size());
dummy = 0;
BasisGridFunction<DefP2Basis,VectorType> def1(Functions::makeFunction(defP2Basis1,dummy));
p2couplings[i] = new LeafP2MortarCoupling<field_type,DeformedGridType>(leafNonmortar, leafMortar);
p2couplings[i]->setup(def0,def1);
}
}
// ///////////////////////////////////////////////////////
// Set up and solve the defect obstacle problem
// ///////////////////////////////////////////////////////
DefectProblem* defectSolver = new DefectProblem(numGrids(), numCouplings());
defectSolver->hessians_ = &matrix;
defectSolver->x_ = &predXNodal_;
defectSolver->rhs_ = &residual;
defectSolver->dirichletNodes_ = &p2DirichletNodes;
// Assemblers for possible rigid obstacles
std::vector<P2OneBodyAssembler<DeformedGridType,field_type>* > rigidObstacles(numCouplings());
for (int i=0; i<numCouplings(); i++) {
CouplingPair<DeformedGridType,DeformedGridType,field_type>& coupling = contactAssembler->coupling_[i];
if (coupling.type_ != CouplingPairBase::RIGID) {
defectSolver->nonmortarLagrangeMatrix_[i] = &(p2couplings[i]->nonmortarLagrangeMatrix_);
defectSolver->mortarLagrangeMatrix_[i] = &(p2couplings[i]->mortarLagrangeMatrix_);
defectSolver->couplings_[i][0] = coupling.gridIdx_[0];
defectSolver->couplings_[i][1] = coupling.gridIdx_[1];
} else {
defectSolver->nonmortarLagrangeMatrix_[i] = NULL;
defectSolver->mortarLagrangeMatrix_[i] = NULL;
defectSolver->couplings_[i][0] = coupling.gridIdx_[0];
defectSolver->couplings_[i][1] = coupling.gridIdx_[0]; // Yes! Twice the same grid!
}
// ///////////////////////////////////
// Set the obstacle values
// ///////////////////////////////////
defectSolver->obstacles_[i].resize(predXNodal_[coupling.gridIdx_[0]].size());
for (size_t j=0; j<defectSolver->obstacles_[i].size(); j++)
defectSolver->obstacles_[i][j].clear();
switch (coupling.type_) {
case CouplingPairBase::NONE:
break;
case CouplingPairBase::CONTACT:
for (size_t j=0; j<defectSolver->obstacles_[i].size(); j++)
if (p2couplings[i]->nonmortarLagrangeMatrix_[j].begin() != p2couplings[i]->nonmortarLagrangeMatrix_[j].end())
defectSolver->obstacles_[i][j].upper(0) = p2couplings[i]->weakObstacle_[j];
break;
case CouplingPairBase::GLUE:
for (size_t j=0; j<defectSolver->obstacles_[i].size(); j++)
if (p2couplings[i]->nonmortarLagrangeMatrix_[j].begin() != p2couplings[i]->nonmortarLagrangeMatrix_[j].end())
for (int k=0; k<dim; k++) {
defectSolver->obstacles_[i][j].lower(k) = 0;
defectSolver->obstacles_[i][j].upper(k) = 0;
}
break;
case CouplingPairBase::RIGID:
rigidObstacles[i] = new P2OneBodyAssembler<DeformedGridType,field_type>(*contactAssembler->grids_[coupling.gridIdx_[0]],
coupling.obsFunction_,
*coupling.patch0(),
coupling.obsDistance_);
rigidObstacles[i]->setObstacles();
defectSolver->obstacles_[i] = rigidObstacles[i]->obstacles_;
break;
default:
DUNE_THROW(Dune::NotImplemented, "Coupling type not implemented!");
}
// /////////////////////////////////////////////////////////////////////////
// Add coordinate system changes for CONTACT and STICK_SLIP couplings
// /////////////////////////////////////////////////////////////////////////
if (coupling.type_==CouplingPairBase::CONTACT
|| coupling.type_ == CouplingPairBase::STICK_SLIP) {
// Make first canonical basis vector
Dune::FieldVector<field_type, dim> e0(0);
e0[0] = 1;
for (size_t j=0; j<defectSolver->obstacles_[i].size(); j++) {
// Modify nonmortar - lagrange matrix
typename MatrixType::row_type::iterator cIt = p2couplings[i]->nonmortarLagrangeMatrix_[j].begin();
typename MatrixType::row_type::iterator cEndIt = p2couplings[i]->nonmortarLagrangeMatrix_[j].end();
for (; cIt!=cEndIt; ++cIt) {
// ///////////////////////////////////////////////
// Householder reflection
// ///////////////////////////////////////////////
BlockType newCoordSystem;
Dune::FieldVector<field_type,dim> field_type_WeakNormal;
for (int k=0; k<dim; k++)
field_type_WeakNormal[k] = p2couplings[i]->weakNormal_[j][k];
householderReflection(e0, field_type_WeakNormal, newCoordSystem);
(*cIt).rightmultiply(newCoordSystem);
}
// Modify nonmortar - lagrange matrix
cIt = p2couplings[i]->mortarLagrangeMatrix_[j].begin();
cEndIt = p2couplings[i]->mortarLagrangeMatrix_[j].end();
for (; cIt!=cEndIt; ++cIt) {
// ///////////////////////////////////////////////
// Householder reflection
// ///////////////////////////////////////////////
BlockType newCoordSystem;
Dune::FieldVector<field_type,dim> field_type_WeakNormal;
for (int k=0; k<dim; k++)
field_type_WeakNormal[k] = p2couplings[i]->weakNormal_[j][k];
householderReflection(e0, field_type_WeakNormal, newCoordSystem);
(*cIt).rightmultiply(newCoordSystem);
}
}
} else if (coupling.type_ == CouplingPairBase::RIGID) {
rigidObstacles[i]->rotateMatrix(matrix[coupling.gridIdx_[0]]);
}
}
// Create a new instance of IpoptApplication
Ipopt::SmartPtr<Ipopt::IpoptApplication> app = new Ipopt::IpoptApplication();
// Change some options
app->Options()->SetNumericValue("tol", 1e-8);
app->Options()->SetNumericValue("acceptable_tol", 1e-8);
app->Options()->SetIntegerValue("max_iter", 100);
app->Options()->SetStringValue("mu_strategy", "adaptive");
app->Options()->SetStringValue("output_file", "ipopt.out");
app->Options()->SetIntegerValue("print_level", 0);
// Specify solver
app->Options()->SetStringValue("linear_solver", "ma27");
// Change memory settings
app->Options()->SetNumericValue("ma27_liw_init_factor", 5);
app->Options()->SetNumericValue("ma27_la_init_factor", 5);
app->Options()->SetNumericValue("ma27_meminc_factor", 10);
// Intialize the IpoptApplication and process the options
Ipopt::ApplicationReturnStatus status;
status = app->Initialize();
if (status != Ipopt::Solve_Succeeded)
DUNE_THROW(SolverError, "Error during IPOpt initialization!");
// Ask Ipopt to solve the problem
Ipopt::SmartPtr<Ipopt::TNLP> defectSolverSmart = defectSolver;
status = app->OptimizeTNLP(defectSolverSmart);
if (status != Ipopt::Solve_Succeeded && status != Ipopt::Solved_To_Acceptable_Level)
DUNE_THROW(SolverError, "Solving the defect problem failed!");
// /////////////////////////////////////////////////////
// Postprocess the solution
// /////////////////////////////////////////////////////
// Transform to canonical coordinates, if necessary
for (int i=0; i<numCouplings(); i++)
if (contactAssembler->coupling_[i].type_ == CouplingPairBase::RIGID)
rigidObstacles[i]->rotateVector(predXNodal_[contactAssembler->coupling_[i].gridIdx_[0]]);
// /////////////////////////////////////////////////
// Change from nodal to hierarchical basis
// /////////////////////////////////////////////////
std::vector<VectorType> hierarchicError(numGrids());
for (int i=0; i<numGrids(); i++) {
VectorType tmp;
P2HierarchicalBasis<typename GridType::LeafGridView, field_type> p2HierarchicalBasis(grids_[i]->leafView());
P2BasisGridFunction predXNodalFunction(*p2NodalBasis[i], predXNodal_[i]);
Functions::interpolate(p2HierarchicalBasis, tmp, predXNodalFunction);
// /////////////////////////////////////////////////////////////////////
// Multiply defect solution in hierarchical basis by the diagonal
// matrix. The edge bubbles can then be used as error indicators.
// /////////////////////////////////////////////////////////////////////
hierarchicError[i].resize(tmp.size());
hierarchicError[i] = 0;
matrix[i].umv(tmp, hierarchicError[i]);
}
// ////////////////////////////////////////////////////////////////
// Set up refinement indicators. Only the edge contributions
// of the hierarchical error are refinement indicators.
// ////////////////////////////////////////////////////////////////
for (int i=0; i<numGrids(); i++) {
P2Basis p2Basis(grids_[i]->leafView());
refinementIndicator[i]->clear();
LeafElementIterator eIt = grids_[i]->template leafbegin<0>();
LeafElementIterator eEndIt = grids_[i]->template leafend<0>();
// Extract all contributions that correspond to edges
for (; eIt!=eEndIt; ++eIt) {
int numLocalDofs = p2Basis.getLocalFiniteElement(*eIt).localCoefficients().size();
for (int j=0; j<numLocalDofs; j++) {
int codim = p2Basis.getLocalFiniteElement(*eIt).localCoefficients().localKey(j).codim();
int subEntity = p2Basis.getLocalFiniteElement(*eIt).localCoefficients().localKey(j).subEntity();
if (codim == dim-1) {
unsigned int globalIdx = p2Basis.index(*eIt, j);
refinementIndicator[i]->set(*eIt, dim-1, subEntity, hierarchicError[i][globalIdx].two_norm());
}
}
}
}
// /////////////////////////////////////////////////
// Clean up
// /////////////////////////////////////////////////
for (int i=0; i<numGrids(); i++) {
delete p2NodalBasis[i];
}
for (int i=0; i<numCouplings(); i++)
delete(p2couplings[i]);
}
template <class GridType, class field_type>
bool HierarchicContactEstimator<GridType,field_type>::DefectProblem::
......
......@@ -373,6 +373,16 @@ public:
std::vector<RefinementIndicator<GridType>* >& refinementIndicator,
std::vector<LocalOperatorAssembler<GridType, P2Lfe, P2Lfe, BlockType>* >& localStiffness);
/** \brief Estimate error for the general large deformation n-body problem */
template <class ProblemType>
void estimate(const std::vector<VectorType>& x,
const std::vector<P1BasisGridFunction*>& volumeTerm,
const std::vector<P1BasisGridFunction*>& neumannTerm,
const std::vector<LeafBoundaryPatch>& dirichletBoundary,
const std::vector<LeafBoundaryPatch>* neumannBoundary,
std::vector<RefinementIndicator<GridType>* >& refinementIndicator,
ProblemType& contactProblem);
/** \brief Compute the error in the energy norm implied by the block diagonal preconditioner. */
field_type getErrorSquared(int grid,
LocalOperatorAssembler<GridType, P2Lfe, P2Lfe, BlockType>* localStiffness) const
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment