From e8570c8e6ab699a64b7d4bbe63ea084bf2a5f9bf Mon Sep 17 00:00:00 2001
From: Jonathan Youett <youett@math.fu-berlin.de>
Date: Mon, 12 Nov 2018 10:44:50 +0100
Subject: [PATCH] One body large deformation contact solver, might need major
cleanup
---
1bnbnonlincontact.cc | 1017 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 1017 insertions(+)
create mode 100644 1bnbnonlincontact.cc
diff --git a/1bnbnonlincontact.cc b/1bnbnonlincontact.cc
new file mode 100644
index 00000000..7b3af24f
--- /dev/null
+++ b/1bnbnonlincontact.cc
@@ -0,0 +1,1017 @@
+#include <config.h>
+
+#include <sys/stat.h>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/sgrid.hh>
+#include <dune/grid/geometrygrid/grid.hh>
+#include <dune/grid/io/file/amirameshwriter.hh>
+#include <dune/grid/io/file/amirameshreader.hh>
+
+#include <dune/istl/io.hh>
+
+#include <dune/fufem/sampleonbitfield.hh>
+#include <dune/fufem/boundarypatchprolongator.hh>
+#include <dune/fufem/readbitfield.hh>
+#include <dune/fufem/estimators/fractionalmarking.hh>
+#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
+#include <dune/fufem/assemblers/operatorassembler.hh>
+#include <dune/fufem/functiontools/gridfunctionadaptor.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+
+#include <dune/solvers/iterationsteps/mmgstep.hh>
+
+#ifdef HAVE_IPOPT
+#include <dune/solvers/solvers/quadraticipopt.hh>
+#endif
+#include <dune/solvers/iterationsteps/trustregiongsstep.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+
+
+// from the 'dune-contact module
+#include <dune/contact/estimators/hierarchiccontactestimator.hh>
+#include <dune/contact/assemblers/onebodyassembler.hh>
+#include <dune/contact/solvers/trustregionstep.hh>
+#include <dune/contact/solvers/contacttransfer.hh>
+#include <dune/contact/solvers/contactobsrestrict.hh>
+#include <dune/contact/solvers/nonlinearcontacttnnmgstep.hh>
+#include <dune/contact/common/onebodystaticcontactproblem.hh>
+
+#include <dune/elasticity/materials/neohookeanmaterial.hh>
+#include <dune/elasticity/assemblers/neohookeassembler.hh>
+#include <dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
+
+#include "../dune-biomech/dune/biomech/amirameshfilehandler.hh"
+
+// The grid dimension
+const int dim = 3;
+
+typedef double field_type;
+
+using namespace Dune;
+using std::string;
+using std::vector;
+using std::map;
+
+// Some types that I need
+typedef BCRSMatrix<FieldMatrix<double, dim, dim> > MatrixType;
+typedef BlockVector<FieldVector<double, dim> > VectorType;
+
+
+double obs(const FieldVector<double, dim>& x) {return 1.00667 + x[dim-1];}
+
+
+class ConstantFunc: public Dune::VirtualFunction<FieldVector<double,dim>,FieldVector<double,dim> >
+{
+public:
+
+ ConstantFunc(FieldVector<double,dim> vec) :
+ vec_(vec)
+ {}
+
+ virtual void evaluate(const FieldVector<double,dim>& x, FieldVector<double,dim>& y) const
+ {
+ y = vec_;
+
+ }
+
+private:
+ FieldVector<double,dim> vec_;
+};
+
+
+
+
+// Concatenate a string and a number
+string makeName(const std::string& name, int number)
+{
+ std::stringstream numberAsAscii;
+ numberAsAscii << number;
+ return name + numberAsAscii.str();
+}
+
+/** \todo The need for the method really is a symptom of bad design */
+size_t index(const std::string& name, const vector<string>& gridNames)
+{
+ for (size_t i=0; i<gridNames.size(); i++)
+ if (name == gridNames[i])
+ return i;
+
+ DUNE_THROW(RangeError, "Grid '" << name << "' not found!");
+}
+
+
+
+int main (int argc, char *argv[]) try
+{
+ // parse data file
+ ParameterTree parameterSet;
+ if (argc==2)
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
+ else
+ ParameterTreeParser::readINITree("1bnbnonlincontact.parset", parameterSet);
+
+ // Some problem settings
+ const int minLevel = parameterSet.get<int>("minLevel");
+ const int maxLevel = parameterSet.get<int>("maxLevel");
+ const bool paramBoundaries = parameterSet.get<bool>("paramBoundaries");
+
+ // Read solver settings
+
+ // _Both_ multigrid solvers
+ const int multigridIterations = parameterSet.get<int>("mgIterations");
+ const int nu1 = parameterSet.get<int>("nu1");
+ const int nu2 = parameterSet.get<int>("nu2");
+ const int mu = parameterSet.get<int>("mu");
+ const int baseIterations = parameterSet.get<int>("baseIt");
+ const double mgTolerance = parameterSet.get<double>("mgTolerance");
+ const double baseTolerance = parameterSet.get<double>("baseTolerance");
+ const double refinementFraction = parameterSet.get("refinementFraction", double(1));
+
+ const bool writeStresses = parameterSet.get<bool>("writeStresses");
+
+ // Trust region parameter
+ const int trustRegionNumIt = parameterSet.get<int>("trustRegionIterations");
+ const double trustRegionTol = parameterSet.get<double>("trustRegionTolerance");
+ double trustRegionRadius = parameterSet.get<double>("trustRegionRadius");
+ const double maxTRRadius = parameterSet.get<double>("maxTRRadius");
+
+ /////////////////////////////////////////////////////////////////////////////////////
+ // Read the actual problem configuration, i.e., grids, and how they are coupled
+ /////////////////////////////////////////////////////////////////////////////////////
+
+ ParameterTree problemConfiguration;
+ ParameterTreeParser::readINITree(parameterSet.get<string>("configurationFile"), problemConfiguration);
+
+ // simulate gravity
+ const bool simulateGravity = problemConfiguration.get<bool>("simulateGravity");
+ const bool useElasticityNorm = problemConfiguration.get<bool>("useLinearElasticityNorm");
+
+ // in which coordinate direction does gravity act
+ const int heightCoord = problemConfiguration.get<int>("heightCoord");
+ // scaling for the Dirichlet values
+ double scaling = parameterSet.get<double>("scaling");
+
+ // read bone-related problem settings
+ string path = problemConfiguration.get<string>("path");
+ string resultPath = problemConfiguration.get<string>("resultPath");
+
+ // Get names of the grids
+ int numGrids = problemConfiguration.get<int>("numGrids");
+ std::vector<string> gridName(numGrids);
+ for (int i=0; i<numGrids; i++) {
+ gridName[i] = problemConfiguration.get<string>(makeName("gridName", i));
+ }
+
+ string gridFile[numGrids];
+ string parFile[numGrids];
+ string dnFile[numGrids];
+ string dvFile[numGrids];
+
+ for (int i=0; i<numGrids; i++) {
+ if (not problemConfiguration.hasSub(gridName[i]))
+ DUNE_THROW(Exception, "No data for '" << gridName[i] << "' is found in the configuration file!");
+ const ParameterTree& gridParameters = problemConfiguration.sub(gridName[i]);
+ gridFile[i] = gridParameters.get<string>("gridFile");
+ parFile[i] = gridParameters.get<string>("parFile");
+ dnFile[i] = gridParameters.get<string>("dnFile");
+ dvFile[i] = gridParameters.get<string>("dvFile");
+ }
+
+ // Read material properties
+ map<string, double> E;
+ map<string, double> Nu;
+ map<string, double> density;
+ map<string, double> D;
+
+ // Read general properties(which will be used if not specified individually)
+ const double gridE_default = problemConfiguration.get("gridE", 1e5);
+ const double gridNu_default = problemConfiguration.get("gridNu", 0.3);
+ const double gridDensity_default = problemConfiguration.get("gridDensity", 1);
+ const double gridD_default = problemConfiguration.get("gridD", 0);
+
+ for (int i=0; i<numGrids; i++) {
+ const ParameterTree& gridParameters = problemConfiguration.sub(gridName[i]);
+ E[gridName[i]] = gridParameters.get("E", gridE_default);
+ Nu[gridName[i]] = gridParameters.get("Nu", gridNu_default);
+ density[gridName[i]] = gridParameters.get("density", gridDensity_default);
+ D[gridName[i]] = gridParameters.get("D", gridD_default);
+ }
+
+ // ////////////////////////////////////////////
+ // Create the various bone grids
+ // ////////////////////////////////////////////
+
+
+ typedef UGGrid<dim> GridType;
+ typedef DeformationFunction<GridType::LeafGridView,VectorType> DeformationFunction;
+ typedef GeometryGrid<GridType,DeformationFunction> DeformedGridType;
+
+
+ typedef BoundaryPatch<GridType::LevelGridView> LevelBoundaryPatch;
+ typedef BoundaryPatch<GridType::LeafGridView> LeafBoundaryPatch;
+
+ map<string,DeformedGridType*> deformedGrids;
+ map<string,GridType*> grids;
+
+ for (int i=0; i<numGrids; i++) {
+
+ grids[gridName[i]] = new GridType;
+
+ if (paramBoundaries)
+ grids[gridName[i]] = AmiraMeshReader<GridType>::read(path + gridFile[i], path + parFile[i]);
+ else
+ grids[gridName[i]] = AmiraMeshReader<GridType>::read(path + gridFile[i]);
+ }
+
+ // initialize the deformed grids with zero displacements
+ std::vector<VectorType> displace(numGrids);
+ std::vector<DeformationFunction*> deformation(numGrids);
+ for (int i=0;i<numGrids;i++) {
+ displace[i].resize(grids[gridName[i]]->size(dim));
+ displace[i]=0;
+ deformation[i]=new DeformationFunction(grids[gridName[i]]->leafView(),displace[i]);
+ deformedGrids[gridName[i]]= new DeformedGridType(grids[gridName[i]],deformation[i]);
+ }
+
+ for (int i=0;i<numGrids;i++)
+ grids[gridName[i]]->setRefinementType(GridType::COPY);
+
+ // ////////////////////////////////////////////////////////////////////
+ // Initial uniform refinement, if requested
+ // ////////////////////////////////////////////////////////////////////
+
+ for (int j=0;j<numGrids;j++) {
+
+ const ParameterTree& gridParameters = problemConfiguration.sub(gridName[j]);
+ // allow different refinement levels for each grid
+ int nRefine = gridParameters.get("nRefines",minLevel);
+
+ for (int i=0; i<nRefine; i++)
+ grids[gridName[j]]->globalRefine(1);
+ }
+
+ // adapt deformed grids
+ for (int i=0;i<numGrids;i++) {
+ displace[i].resize(grids[gridName[i]]->size(dim));
+ displace[i]=0;
+ deformation[i]->setDeformation(displace[i]);
+ deformedGrids[gridName[i]]->update();
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+ // Set up hierarchy of bitfields containing the Dirichlet flags
+ //////////////////////////////////////////////////////////////////////////
+ vector<VectorType> coarseDirichletValues(numGrids);
+
+ for (int i=0; i<numGrids; i++) {
+
+ coarseDirichletValues[i].resize(grids[gridName[i]]->size(0, dim));
+ coarseDirichletValues[i] = 0;
+
+ AmiraMeshReader<GridType>::readFunction(coarseDirichletValues[i], path + dvFile[i]);
+
+ // Scale Dirichlet values
+ coarseDirichletValues[i] *= scaling;
+ }
+
+ // Set up hierarchy of bitfields containing the Dirichlet flags
+ vector<LevelBoundaryPatch> coarseDirichletBoundary(numGrids);
+ vector<LeafBoundaryPatch> leafDirichletBoundary(numGrids);
+
+ for (int i=0; i<numGrids; i++) {
+ coarseDirichletBoundary[i].setup(grids[gridName[i]]->levelView(0));
+ leafDirichletBoundary[i].setup(grids[gridName[i]]->leafView());
+ readBoundaryPatch<GridType>(coarseDirichletBoundary[i], path + dnFile[i]);
+ }
+
+ vector<BitSetVector<dim> > dirichletNodes(numGrids);
+
+ // Get overall top level
+ int toplevel = 0;
+ for (int i=0;i<numGrids;i++)
+ toplevel = std::max(toplevel, grids[gridName[i]]->maxLevel());
+
+ std::vector<VectorType> dirichletValues(numGrids);
+
+ for (int i=0; i<numGrids; i++) {
+
+ // Extend the dirichlet information to all grid levels
+ BoundaryPatchProlongator<GridType>::prolong(coarseDirichletBoundary[i], leafDirichletBoundary[i]);
+
+ // Copy vertices into bitfield
+ int fSSize = grids[gridName[i]]->size(dim);
+ dirichletNodes[i].resize(fSSize, false);
+ for (int k=0; k<fSSize; k++)
+ dirichletNodes[i][k] = leafDirichletBoundary[i].containsVertex(k);
+
+ sampleOnBitField(*grids[gridName[i]], coarseDirichletValues[i], dirichletValues[i], dirichletNodes[i]);
+ std::cout<<dirichletNodes[i].count()<<" Dirichlet nodes for grid "<<i<<std::endl;
+ }
+
+ // displacements of continua
+ vector<VectorType> u(numGrids),dirVals(numGrids), extForces(numGrids);
+
+ // add Dirichlet values
+ for (int i=0; i<numGrids; i++) {
+ u[i].resize(grids[gridName[i]]->size(dim));
+ u[i] = 0;
+ extForces[i] = u[i];
+ dirVals[i] = u[i];
+ for (size_t j=0; j<dirVals[i].size(); j++)
+ for (int k=0; k<dim; k++)
+ if (dirichletNodes[i][j][k])
+ dirVals[i][j][k] = dirichletValues[i][j][k];
+ }
+
+ if (problemConfiguration.get<bool>("restart")) {
+ int loadStep = problemConfiguration.get<int>("loadingStep");
+
+ for (int i=0; i<numGrids; i++) {
+
+ std::string restartFile = genFilename(resultPath, makeName("grid",i) + "loading", loadStep);
+ AmiraMeshReader<GridType>::readFunction(u[i], restartFile);
+
+ }
+ }
+
+
+ // make dirichlet bitfields containing dirichlet information for all grids
+ int size = 0;
+ int offset = 0;
+
+ for (int i=0; i<numGrids; i++)
+ size += dirichletValues[i].size();
+
+ BitSetVector<dim> totalDirichletNodes(size);
+
+ for (int i=0; i<numGrids; i++) {
+ for (size_t j=0; j<dirichletValues[i].size(); j++)
+ totalDirichletNodes[offset + j] = dirichletNodes[i][j];
+ offset += dirichletValues[i].size();
+ }
+
+ // ///////////////////////////////////
+ // Create contact assembler
+ // ///////////////////////////////////
+
+ // Set mortar couplings. Currently not more than one for each grid
+ BoundaryPatch<DeformedGridType::LevelGridView> nonmortarPatch;
+ const ParameterTree& couplingParameters = problemConfiguration.sub(makeName("coupling",0));
+
+ // the involved grids
+ string nonMortarGrid = couplingParameters.get<string>("nonmortarGrid");
+
+ // read field describing the nonmortar boundary
+ nonmortarPatch.setup(deformedGrids[gridName[0]]->levelView(0),true);
+ readBoundaryPatch<DeformedGridType>(nonmortarPatch, path + couplingParameters.get<string>("nonmortarPatch"));
+
+ OneBodyAssembler<DeformedGridType,VectorType> contactAssembler(*deformedGrids[gridName[0]], obs, nonmortarPatch, couplingParameters.get<double>("couplingDistance"));
+
+
+ std::vector<BitSetVector<1> > allHasObs(toplevel+1);
+ allHasObs.back().resize(size,true);
+
+ for (size_t i=0; i<totalDirichletNodes.size(); i++)
+ if (totalDirichletNodes[i][0])
+ allHasObs.back()[i].flip();
+
+ for (size_t i=0; i<grids[gridName[0]]->maxLevel()+1; i++)
+ allHasObs[i].resize(grids[gridName[0]]->size(i,dim),true);
+
+ // ////////////////////////////////////////////////////
+ // Create a multigrid solver and trust region solver
+ // ////////////////////////////////////////////////////
+
+ EnergyNorm<MatrixType, VectorType> elastNorm;
+
+ // First create a base solver
+#ifndef HAVE_IPOPT
+#error "You need IpOpt!"
+#endif
+ /*
+ TrustRegionGSStep<MatrixType, VectorType> baseSolverStep;
+
+ ::LoopSolver<VectorType> baseSolver(&baseSolverStep,
+ baseIterations,
+ baseTolerance,
+ &elastNorm,
+ Solver::REDUCED);
+ */
+
+
+ QuadraticIPOptSolver<MatrixType,VectorType> baseSolver;
+ baseSolver.verbosity_ = Solver::QUIET;
+ baseSolver.tolerance_ = baseTolerance;
+
+
+
+ // Make pre and postsmoothers
+ TrustRegionGSStep<MatrixType, VectorType> presmoother, postsmoother;
+
+ MonotoneMGStep<MatrixType, VectorType> multigridStep;
+
+ multigridStep.setMGType(mu, nu1, nu2);
+ multigridStep.basesolver_ = &baseSolver;
+ multigridStep.setSmoother(&presmoother, &postsmoother);
+ multigridStep.obstacleRestrictor_ = new ContactObsRestriction<VectorType>;
+ multigridStep.hasObstacle_ = &allHasObs;
+
+ multigridStep.ignoreNodes_ = &totalDirichletNodes;
+
+ // transfer operator need to be set every Newton iteration
+ std::vector<ContactMGTransfer<VectorType>* > mgTransfers(toplevel);
+ for (size_t i=0; i<mgTransfers.size(); i++)
+ mgTransfers[i] = new ContactMGTransfer<VectorType>;
+
+ multigridStep.setTransferOperators(mgTransfers);
+
+ EnergyNorm<MatrixType, VectorType > energyNorm(multigridStep);
+
+ ::LoopSolver<VectorType> solver(&multigridStep,
+ multigridIterations,
+ mgTolerance,
+ (useElasticityNorm) ? &elastNorm : &energyNorm,
+ Solver::FULL);
+
+ // Create the materials
+ typedef P1NodalBasis<GridType::LeafGridView> P1Basis;
+ std::vector<P1Basis*> coarseP1Bases(numGrids);
+ for (int j=0; j<numGrids; j++)
+ coarseP1Bases[j] = new P1Basis(grids[gridName[j]]->leafView());
+
+ typedef NeoHookeMaterial<P1Basis> MaterialType;
+ //typedef GeomExactStVenantMaterial<P1Basis> MaterialType;
+ std::vector<MaterialType> materials(numGrids);
+ for (int i=0; i<numGrids;i++)
+ materials[i].setup(E[gridName[i]],Nu[gridName[i]],*coarseP1Bases[i]);
+
+ // Make static contact problem
+ typedef OneBodyStaticContactProblem<VectorType,MatrixType, MaterialType> ContactProblem;
+ ContactProblem contactProblem(contactAssembler, materials, extForces);
+
+ // Create trust-region step
+ TrustRegionStep<ContactProblem> trustRegionStep(contactProblem, solver, trustRegionRadius, maxTRRadius);
+
+ // ///////////////////////////////////////////////////
+ // Do a homotopy of the Dirichlet boundary data
+ // ///////////////////////////////////////////////////
+
+ // //////////////////////////////////////////////////////////
+ // Assemble linear elasticity matrix for the energy norm
+ // used by the MMG termination criterion
+ // //////////////////////////////////////////////////////////
+
+ std::vector<MatrixType*> elast(numGrids);
+ MatrixType bigElastMatrix;
+
+ if (useElasticityNorm) {
+ for (int i=0; i<numGrids; i++) {
+
+ OperatorAssembler<P1Basis,P1Basis> globalAssembler(*coarseP1Bases[i],*coarseP1Bases[i]);
+
+ StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> localAssembler(E[gridName[i]],Nu[gridName[i]]);
+
+ elast[i] = new MatrixType;
+ globalAssembler.assemble(localAssembler, *elast[i]);
+ }
+
+
+ MatrixIndexSet matInd(size,size);
+ size_t matOffSet = 0;
+ for (int i=0; i<numGrids; i++) {
+
+ for (size_t iRow = 0; iRow<elast[i]->N(); iRow++) {
+
+ const MatrixType::row_type& row = (*elast[i])[iRow];
+
+ // Loop over all columns of the stiffness matrix
+ MatrixType::row_type::const_iterator j = row.begin();
+ MatrixType::row_type::const_iterator jEnd = row.end();
+
+ for (; j!=jEnd; ++j)
+ matInd.add(iRow+matOffSet, matOffSet+j.index());
+ }
+
+ matOffSet += u[i].size();
+ }
+
+ matInd.exportIdx(bigElastMatrix);
+ bigElastMatrix = 0;
+ matOffSet = 0;
+
+ for (int i=0; i<numGrids; i++) {
+
+ for (size_t iRow = 0; iRow<elast[i]->N(); iRow++) {
+
+ const MatrixType::row_type& row = (*elast[i])[iRow];
+
+ // Loop over all columns of the stiffness matrix
+ MatrixType::row_type::const_iterator j = row.begin();
+ MatrixType::row_type::const_iterator jEnd = row.end();
+
+ for (; j!=jEnd; ++j)
+ bigElastMatrix[iRow+matOffSet][matOffSet+j.index()] = *j;
+ }
+
+ matOffSet += u[i].size();
+ }
+
+ elastNorm.setMatrix(&bigElastMatrix);
+ }
+
+
+ double loadFactor = 0;
+ double loadIncrement = problemConfiguration.get<double>("loadIncrement");
+
+ int counter =0;
+
+ if (problemConfiguration.get<bool>("restart")) {
+ counter = problemConfiguration.get<int>("loadingStep") +1;
+ loadFactor = problemConfiguration.get<double>("loadFactor");
+ loadIncrement = problemConfiguration.get<double>("loadIncrement");
+ }
+
+ double avMgIt = 0; double infeasibility = 1;
+ double nonlinInfeasibility = 1;
+ do {
+ /*
+ if (loadFactor < 1) {
+
+ // ////////////////////////////////////////////////////
+ // Compute new feasible load increment
+ // ////////////////////////////////////////////////////
+
+ std::vector<VectorType> uSave=u;
+ loadIncrement *= 2; // double it once, cause we're going to half it at least once, too!
+ bool isNan;
+ do {
+
+ loadIncrement /= 2;
+
+ // Set new Dirichlet values in solution
+ for (int j=0; j<numGrids; j++)
+ for (int k=0; k<u[j].size(); k++)
+ for (int l=0; l<dim; l++)
+ if (dirichletNodes[j][k][l])
+ uSave[j][k][l] = dirichletValues[j][k][l] * (loadFactor + loadIncrement);
+
+ isNan = false;
+ for (int i=0; i<numGrids; i++) {
+ double eng = materials[i].energy(uSave[i]);
+ std::cout<<"energy "<<eng<<std::endl;
+ isNan |= isnan(eng);
+ }
+ } while (isNan);
+
+ loadFactor += loadIncrement;
+
+ std::cout << "####################################################" << std::endl;
+ std::cout << "New load factor: " << loadFactor
+ << " new load increment: " << loadIncrement << std::endl;
+ std::cout << "####################################################" << std::endl;
+
+ } */
+ loadFactor += loadIncrement;
+
+ std::vector<VectorType> uSave=u;
+ // Set new Dirichlet values in solution
+ for (int j=0; j<numGrids; j++)
+ for (int k=0; k<u[j].size(); k++)
+ for (int l=0; l<dim; l++)
+ if (dirichletNodes[j][k][l])
+ uSave[j][k][l] = dirichletValues[j][k][l] * (loadFactor);
+
+ for (int i=0; i<numGrids; i++) {
+ double eng = materials[i].energy(uSave[i]);
+ std::cout<<"energy "<<eng<<std::endl;
+ }
+ std::cout << "####################################################" << std::endl;
+ std::cout << "New load factor: " << loadFactor
+ << " new load increment: " << loadIncrement << std::endl;
+ std::cout << "####################################################" << std::endl;
+
+
+
+ // Set new Dirichlet values in solution
+ for (int j=0; j<numGrids; j++)
+ for (int k=0; k<u[j].size(); k++)
+ for (int l=0; l<dim; l++)
+ if (dirichletNodes[j][k][l])
+ //dirVals[j][k][l] = dirichletValues[j][k][l] * (loadFactor);
+ dirVals[j][k][l] = dirichletValues[j][k][l] * (loadIncrement);
+
+ // /////////////////////////////////////////////////////
+ // Solve coarse problem
+ // /////////////////////////////////////////////////////
+
+ for (int j=0; j<numGrids; j++)
+ for (int k=0; k<u[j].size(); k++)
+ if (dirichletNodes[j][k][0]) {
+ std::cout<<"Dir "<<dirVals[j][k]<<std::endl;
+ break;
+ }
+
+ std::ostringstream s;
+ int zer(0);
+ s << counter;
+
+ std::string itepath = resultPath + "mergedGrids";
+ int new_dir2 = mkdir(itepath.c_str(), 0744);
+ itepath += "/"+s.str()+"/";
+ //Create loging files for several entities
+ // new folder for results
+ new_dir2 = mkdir(itepath.c_str(), 0744);
+ if (new_dir2 != -1)
+ printf("Created new folder for results\n");
+ else
+ printf("Error while creating folder!\n");
+ // Newton iteration for linearized equation
+
+
+ // actual trust region iteration and error
+ int trIt = 0;
+ double error = trustRegionTol+1;
+ double relError = trustRegionTol+1;
+
+ LeafBoundaryPatch allBound(grids[gridName[0]]->leafView(),true);
+ std::vector<int> allGloToLo;
+ allBound.makeGlobalToLocal(allGloToLo);
+
+ while ((relError >= trustRegionTol) && trIt < trustRegionNumIt) {
+ ++trIt;
+
+ // Deform grids with actual displacement iterate
+ for (int i=0;i<numGrids;i++)
+ deformation[i]->setDeformation(u[i]);
+
+
+ // compute obstacles and setup the MMG solver
+ contactAssembler.setObstacles();
+ multigridStep.obstacles_ = &contactAssembler.obstacles_;
+
+ std::vector<int> gl;
+ nonmortarPatch.makeGlobalToLocal(gl);
+
+ std::cout<<"obstacles \n";
+ for (int i=0; i<u[0].size();i++)
+ if (gl[i]>-1)
+ std::cout<<contactAssembler.obstacles_.back()[i].lower(dim-1)<<std::endl;
+
+ for (size_t i=0; i<mgTransfers.size(); i++) {
+ mgTransfers[i]->setup(*deformedGrids[gridName[0]], i, i+1,
+ contactAssembler.localCoordSystems_[i],
+ contactAssembler.localCoordSystems_[i+1]);
+ }
+
+ //transform iterate into mortar basis coordinates
+ VectorType totalIterate = u[0];
+ VectorType newTotalIterate = totalIterate;
+ contactAssembler.rotateVector(newTotalIterate);
+ // setup trust region problem
+ if (trIt==1)
+ trustRegionStep.setProblem(newTotalIterate, contactProblem, dirVals, true);
+ else
+ trustRegionStep.setProblem(newTotalIterate, contactProblem, extForces, false);
+
+ trustRegionStep.iterate();
+ avMgIt += solver.ite_;
+
+ // get new iterate
+ newTotalIterate = trustRegionStep.getSol();
+
+ // transform back into nodal coordinates
+ contactAssembler.rotateVector(newTotalIterate);
+ u[0] = newTotalIterate;
+
+ for (int i=0; i<numGrids; i++) {
+
+ std::string filme2 = genFilename(itepath, makeName("grid",i) + "iteration", trIt);
+
+ LeafAmiraMeshWriter<GridType> amiramesh;
+ amiramesh.addLeafGrid(*grids[gridName[i]],true);
+ amiramesh.addVertexData(u[i], grids[gridName[i]]->leafView());
+ amiramesh.write(filme2);
+ }
+
+
+ // get the correction
+ newTotalIterate -= totalIterate;
+ std::vector<VectorType> correction(numGrids);
+ correction[0] = newTotalIterate;
+
+ // compute error estimate
+ error = 0;
+ double error2 = 0;
+
+ for (int i=0;i<numGrids;i++) {
+ if (useElasticityNorm)
+ error += EnergyNorm<MatrixType,VectorType>::normSquared(correction[i], *elast[i]);
+ else
+ error += EnergyNorm<MatrixType,VectorType>::normSquared(correction[i], *contactProblem.stiffMat_[i]);
+ if (useElasticityNorm)
+ error2 += EnergyNorm<MatrixType,VectorType>::normSquared(u[i], *elast[i]);
+ else
+ error2 += EnergyNorm<MatrixType,VectorType>::normSquared(u[i], *contactProblem.stiffMat_[i]);
+ }
+
+ relError = error/error2;
+
+ std::cout << "In Trust Region iteration "<<trIt<<" , absolut norm of correction : "<<error<<
+ " relative error "<<relError<<"\n";
+
+ }
+
+ avMgIt /= trIt;
+
+ // Output result
+ for (int i=0; i<numGrids; i++) {
+
+ std::string filename = genFilename(resultPath, makeName("grid",i) + "loading", counter);
+ LeafAmiraMeshWriter<GridType>::writeBlockVector(*grids[gridName[i]], u[i], filename,true);
+ }
+
+ counter++;
+
+ } while (loadFactor < 1);
+
+ std::cout<<"--------------------------------------------------- \n";
+ std::cout<<"| Finished Loading phase \n";
+ std::cout<<"--------------------------------------------------- \n";
+
+ std::cout<<"average multigrid iterations are "<<avMgIt<<std::endl;
+
+ // /////////////////////////////////////////////////////////////////////
+ // Refinement Loop
+ // /////////////////////////////////////////////////////////////////////
+/*
+ for (int level=1; level<=maxLevel; level++) {
+
+ // ////////////////////////////////////////////////////////////////////////////
+ // Refine locally and transfer the current solution to the new leaf level
+ // ////////////////////////////////////////////////////////////////////////////
+
+ // GeometricEstimator<GridType> estimator;
+
+ // estimator.estimate(grid[0], (toplevel<=minLevel) ? refineAll : refineCondition, x[0]);
+ // estimator.estimate(grid[1], (toplevel<=minLevel) ? refineAll : refineCondition, x[1]);
+
+ HierarchicContactEstimator<GridType> estimator(grid[0], grid[1]);
+
+ estimator.couplings_[0].obsPatch_ = &untransfObsPatch;
+ estimator.couplings_[0].obsDistance_ = obsDistance;
+ estimator.couplings_[0].type_ = CouplingPairBase::CONTACT;
+
+ array<LeafBoundaryPatch, 2> leafDirichletBoundary;
+ for (int j=0; j<grid.size(); j++)
+ BoundaryPatchProlongator<GridType>::prolong(coarseDirichletBoundary[j], leafDirichletBoundary[j]);
+
+ OgdenMaterialLocalStiffness<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> localOgdenStiffness(E, nu, d);
+
+ std::vector<RefinementIndicator<GridType>*> refinementIndicator(2);
+ refinementIndicator[0] = new RefinementIndicator<GridType>(grid[0]);
+ refinementIndicator[1] = new RefinementIndicator<GridType>(grid[1]);
+
+ estimator.estimate(x[0], x[1],
+ &leafDirichletBoundary[0],
+ &leafDirichletBoundary[1],
+ refinementIndicator,
+ &localOgdenStiffness, &localOgdenStiffness);
+
+ // ////////////////////////////////////////////////////
+ // Refine grids
+ // ////////////////////////////////////////////////////
+
+ std::vector<GridType*> adaptiveGridVector(2);
+ adaptiveGridVector[0] = &grid[0];
+ adaptiveGridVector[1] = &grid[1];
+ FractionalMarkingStrategy<GridType>::mark(refinementIndicator, adaptiveGridVector, refinementFraction);
+
+ for (int i=0; i<2; i++) {
+
+ P1Basis p1Basis(grid[i].leafView());
+ GridFunctionAdaptor<P1Basis> adaptor(p1Basis,true,true);
+
+ grid[i].preAdapt();
+ grid[i].adapt();
+ grid[i].postAdapt();
+
+ p1Basis.update();
+ adaptor.adapt(x[i]);
+
+ deform[i]->updateDeformation(x[i],grid[i].maxLevel());
+ defGrid[i]->update();
+
+ }
+
+ std::cout << "########################################################" << std::endl;
+ std::cout << " Grids refined" << std::endl;
+ std::cout << " Grid: 0 Level: " << grid[0].maxLevel()
+ << " vertices: " << grid[0].size(grid[0].maxLevel(), dim)
+ << " elements: " << grid[0].size(grid[0].maxLevel(), 0) << std::endl;
+ std::cout << " Grid: 1 Level: " << grid[1].maxLevel()
+ << " vertices: " << grid[1].size(grid[1].maxLevel(), dim)
+ << " elements: " << grid[1].size(grid[1].maxLevel(), 0) << std::endl;
+ std::cout << "########################################################" << std::endl;
+
+ std::cout << "####################################################" << std::endl;
+ std::cout << " Solving on level: " << toplevel << std::endl;
+ std::cout << "####################################################" << std::endl;
+
+ // does this keep the level 0 entries?
+ dirichletBoundary[0].resize(toplevel+1);
+ dirichletBoundary[1].resize(toplevel+1);
+ dirichletBoundary[0][0] = coarseDirichletBoundary[0];
+ dirichletBoundary[1][0] = coarseDirichletBoundary[1];
+
+
+ dirichletNodes[0].resize(toplevel+1);
+ dirichletNodes[1].resize(toplevel+1);
+
+ BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary[0]);
+ BoundaryPatchProlongator<GridType>::prolong(dirichletBoundary[1]);
+
+ for (int i=0; i<=toplevel; i++) {
+
+ int fSSize0 = grid[0].size(i,dim);
+ dirichletNodes[0][i].resize(fSSize0*dim);
+ for (int j=0; j<fSSize0; j++)
+ for (int k=0; k<dim; k++)
+ dirichletNodes[0][i][j*dim+k] = dirichletBoundary[0][i].containsVertex(j);
+
+ int fSSize1 = grid[1].size(i,dim);
+ dirichletNodes[1][i].resize(fSSize1 * dim);
+ for (int j=0; j<fSSize1; j++)
+ for (int k=0; k<dim; k++)
+ dirichletNodes[1][i][j*dim+k] = dirichletBoundary[1][i].containsVertex(j);
+
+ }
+
+ dirichletValues[0].resize(toplevel+1);
+ dirichletValues[1].resize(toplevel+1);
+ dirichletValues[0][0] = coarseDirichletValues[0];
+ dirichletValues[1][0] = coarseDirichletValues[1];
+
+ sampleOnBitField(grid[0], dirichletValues[0], dirichletNodes[0]);
+ sampleOnBitField(grid[1], dirichletValues[1], dirichletNodes[1]);
+
+ // ///////////////////////////////////////////////
+ // Create new rhs vectors
+ // ///////////////////////////////////////////////
+
+ rhs[0].resize(grid[0].size(toplevel,dim));
+ rhs[1].resize(grid[1].size(toplevel,dim));
+
+ // Set right hand side vectors
+ for (int i=0; i<2; i++) {
+
+ rhs[i] = 0;
+
+ for (int j=0; j<rhs[i].size(); j++)
+ for (int k=0; k<dim; k++) {
+ if (dirichletNodes[i][toplevel][j][k])
+ x[i][j][k] = dirichletValues[i][toplevel][j][k];
+
+ }
+
+ }
+
+ // /////////////////////////////////////////////////////////////////////////////
+ // Make dirichlet bitfields containing dirichlet information for both grids
+ // /////////////////////////////////////////////////////////////////////////////
+ totalDirichletNodes.resize(toplevel+1);
+
+ for (int i=0; i<=toplevel; i++) {
+
+ int offset = dirichletValues[0][i].size();
+
+ totalDirichletNodes[i].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];
+
+ for (int j=0; j<dirichletNodes[1][i].size(); j++)
+ totalDirichletNodes[i][offset + j] = dirichletNodes[1][i][j];
+
+
+ }
+
+ trustRegionObstacles.resize(toplevel+1);
+ for (int i=0; i<=toplevel; i++)
+ trustRegionObstacles[i].resize(grid[0].size(i, dim) + grid[1].size(i, dim));
+
+
+ // /////////////////////////////////////////////////////
+ // Assemble the obstacle and the transfer operator
+ // /////////////////////////////////////////////////////
+
+ //update the deformed grids!
+ for (int j=0;j<2;j++)
+ deform[j]->updateDeformation(x[j],grid[j].maxLevel());
+
+ contactAssembler.assembleObstacle();
+ contactAssembler.assembleTransferOperator();
+
+ // //////////////////////////////////////////////////////
+ // Lengthen obstacle vector, to formally account for both
+ // bodies, even though the second one doesn't actually have obstacles
+ // //////////////////////////////////////////////////////
+
+ hasObstacle.resize(toplevel+1);
+ for (int j=0; j<=toplevel; j++)
+ hasObstacle[j].resize(grid[0].size(j, dim) + grid[1].size(j,dim), true);
+
+ // //////////////////////////////////////////////////////////////////////
+ // The MG transfer operators involve the mortar coupling
+ // coupling operator and have to be reassembled at each loading step
+ // //////////////////////////////////////////////////////////////////////
+ for (int k=0; k<multigridStep.mgTransfer_.size(); k++)
+ delete(multigridStep.mgTransfer_[k]);
+
+ multigridStep.mgTransfer_.resize(toplevel);
+
+ for (int k=0; k<multigridStep.mgTransfer_.size(); k++) {
+
+ ContactMGTransfer<VectorType>* newTransferOp = new ContactMGTransfer<VectorType>;
+ newTransferOp->setup(grid[0], grid[1], k, k+1,
+ contactAssembler.contactCoupling_.mortarLagrangeMatrix(k),
+ contactAssembler.localCoordSystems_[k],
+ contactAssembler.localCoordSystems_[k+1],
+ *contactAssembler.nonmortarBoundary_[0][k].getVertices(),
+ *contactAssembler.nonmortarBoundary_[0][k+1].getVertices());
+
+ multigridStep.mgTransfer_[k] = newTransferOp;
+
+
+ }
+
+ // //////////////////////////////////////////////////////////
+ // Assemble linear elasticity matrix for the energy norm
+ // used by the MMG termination criterion
+ // //////////////////////////////////////////////////////////
+ // Assemble separate linear elasticity problems
+ P1Basis p1Basis0(grid[0].leafView());
+ P1Basis p1Basis1(grid[1].leafView());
+
+ OperatorAssembler<P1Basis,P1Basis> globalAssembler0(p1Basis0,p1Basis0);
+ OperatorAssembler<P1Basis,P1Basis> globalAssembler1(p1Basis1,p1Basis1);
+
+ StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement, P1Basis::LocalFiniteElement> localAssembler(2.5e5, 0.3);
+
+ MatrixType stiffnessMatrix0, stiffnessMatrix1;
+ globalAssembler0.assemble(localAssembler, stiffnessMatrix0);
+ globalAssembler1.assemble(localAssembler, stiffnessMatrix1);
+
+ array<const MatrixType*, 2> submat;
+ submat[0] = &stiffnessMatrix0;
+ submat[1] = &stiffnessMatrix1;
+
+ MatrixType transformedLinearMatrix;
+ contactAssembler.assemble(submat, transformedLinearMatrix);
+
+ dynamic_cast<EnergyNorm<MatrixType,VectorType>*>(solver.errorNorm_)->setMatrix(&transformedLinearMatrix);
+
+ // /////////////////////////////////////////////////////
+ // Trust-Region Solver
+ // /////////////////////////////////////////////////////
+
+ solveTrustRegionProblem<DefGridType, MatrixType, VectorType, P1Basis, P1Basis>(ogdenAssembler,
+ contactAssembler,
+ x, rhs, grid[0].maxLevel(),
+ trustRegionRadius,
+ trustRegionObstacles,
+ totalDirichletNodes,
+ solver,
+ multigridStep,
+ maxNewtonStepsII);
+
+
+ // ///////////////////////
+ // Output result
+ // ///////////////////////
+ LeafAmiraMeshWriter<GridType> amiramesh0;
+ amiramesh0.addLeafGrid(grid[0],true);
+ amiramesh0.addVertexData(x[0], grid[0].leafView());
+ amiramesh0.write("0resultGrid");
+
+ LeafAmiraMeshWriter<GridType> amiramesh1;
+ amiramesh1.addLeafGrid(grid[1],true);
+ amiramesh1.addVertexData(x[1], grid[1].leafView());
+ amiramesh1.write("1resultGrid");
+
+ }
+ */
+} catch (Exception e) {
+
+ std::cout << e << std::endl;
+
+}
+
+
+
+
--
GitLab