From 2dd7af0200d288e026d23516eba467c8a6738ff0 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 23 Mar 2015 08:01:46 +0000
Subject: [PATCH] Make compile with the new dune-functions FE bases
[[Imported from SVN: r10102]]
---
src/mixed-cosserat-continuum.cc | 64 ++++++++++++++++-------------
src/rod-eoc.cc | 16 +++++++-
src/rodobstacle.cc | 72 +++++++++++++++++----------------
3 files changed, 87 insertions(+), 65 deletions(-)
diff --git a/src/mixed-cosserat-continuum.cc b/src/mixed-cosserat-continuum.cc
index bce35209..d96af1ef 100644
--- a/src/mixed-cosserat-continuum.cc
+++ b/src/mixed-cosserat-continuum.cc
@@ -25,11 +25,12 @@
#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>
#include <dune/fufem/functiontools/basisinterpolator.hh>
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/functionspacebases/p2nodalbasis.hh>
+#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
#include <dune/fufem/dunepython.hh>
#include <dune/solvers/solvers/iterativesolver.hh>
@@ -246,13 +247,22 @@ int main (int argc, char *argv[]) try
typedef GridType::LeafGridView GridView;
GridView gridView = grid->leafGridView();
- typedef P3NodalBasis<GridView,double> DeformationFEBasis;
- typedef P2NodalBasis<GridView,double> OrientationFEBasis;
+ typedef Dune::Functions::PQKNodalBasis<GridView,3> DeformationFEBasis;
+ typedef Dune::Functions::PQKNodalBasis<GridView,2> OrientationFEBasis;
DeformationFEBasis deformationFEBasis(gridView);
OrientationFEBasis orientationFEBasis(gridView);
- std::cout << "Deformation: " << deformationFEBasis.size() << ", orientation: " << orientationFEBasis.size() << std::endl;
+ // Construct fufem-style function space bases to ease the transition to dune-functions
+ typedef DuneFunctionsBasis<DeformationFEBasis> FufemDeformationFEBasis;
+ FufemDeformationFEBasis fufemDeformationFEBasis(deformationFEBasis);
+
+ typedef DuneFunctionsBasis<OrientationFEBasis> FufemOrientationFEBasis;
+ FufemOrientationFEBasis fufemOrientationFEBasis(orientationFEBasis);
+
+
+
+ std::cout << "Deformation: " << deformationFEBasis.indexSet().size() << ", orientation: " << orientationFEBasis.indexSet().size() << std::endl;
// /////////////////////////////////////////
// Read Dirichlet values
@@ -294,20 +304,20 @@ int main (int argc, char *argv[]) try
std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
- BitSetVector<1> deformationDirichletNodes(deformationFEBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,deformationFEBasis,deformationDirichletNodes);
+ BitSetVector<1> deformationDirichletNodes(deformationFEBasis.indexSet().size(), false);
+ constructBoundaryDofs(dirichletBoundary,fufemDeformationFEBasis,deformationDirichletNodes);
- BitSetVector<3> deformationDirichletDofs(deformationFEBasis.size(), false);
- for (size_t i=0; i<deformationFEBasis.size(); i++)
+ BitSetVector<3> deformationDirichletDofs(deformationFEBasis.indexSet().size(), false);
+ for (size_t i=0; i<deformationFEBasis.indexSet().size(); i++)
if (deformationDirichletNodes[i][0])
for (int j=0; j<3; j++)
deformationDirichletDofs[i][j] = true;
- BitSetVector<1> orientationDirichletNodes(orientationFEBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,orientationFEBasis,orientationDirichletNodes);
+ BitSetVector<1> orientationDirichletNodes(orientationFEBasis.indexSet().size(), false);
+ constructBoundaryDofs(dirichletBoundary,fufemOrientationFEBasis,orientationDirichletNodes);
- BitSetVector<3> orientationDirichletDofs(orientationFEBasis.size(), false);
- for (size_t i=0; i<orientationFEBasis.size(); i++)
+ BitSetVector<3> orientationDirichletDofs(orientationFEBasis.indexSet().size(), false);
+ for (size_t i=0; i<orientationFEBasis.indexSet().size(); i++)
if (orientationDirichletNodes[i][0])
for (int j=0; j<3; j++)
orientationDirichletDofs[i][j] = true;
@@ -316,18 +326,18 @@ int main (int argc, char *argv[]) try
// Initial iterate
// //////////////////////////
- DeformationSolutionType xDisp(deformationFEBasis.size());
+ DeformationSolutionType xDisp(deformationFEBasis.indexSet().size());
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
std::vector<FieldVector<double,3> > v;
- Functions::interpolate(deformationFEBasis, v, pythonInitialDeformation);
+ ::Functions::interpolate(fufemDeformationFEBasis, v, pythonInitialDeformation);
for (size_t i=0; i<xDisp.size(); i++)
xDisp[i] = v[i];
- OrientationSolutionType xOrient(orientationFEBasis.size());
+ OrientationSolutionType xOrient(orientationFEBasis.indexSet().size());
#if 0
lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
@@ -343,8 +353,8 @@ int main (int argc, char *argv[]) try
////////////////////////////////////////////////////////
// Output initial iterate (of homotopy loop)
- CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,xDisp,
- orientationFEBasis,xOrient,
+ CosseratVTKWriter<GridType>::writeMixed<FufemDeformationFEBasis,FufemOrientationFEBasis>(fufemDeformationFEBasis,xDisp,
+ fufemOrientationFEBasis,xOrient,
resultPath + "mixed-cosserat_homotopy_0");
for (int i=0; i<numHomotopySteps; i++) {
@@ -372,17 +382,15 @@ int main (int argc, char *argv[]) try
}
// Assembler using ADOL-C
- MixedCosseratEnergy<GridView,
- DeformationFEBasis::LocalFiniteElement,
- OrientationFEBasis::LocalFiniteElement,
+ MixedCosseratEnergy<DeformationFEBasis,
+ OrientationFEBasis,
3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters,
&neumannBoundary,
neumannFunction.get());
- MixedLocalGFEADOLCStiffness<GridView,
- DeformationFEBasis::LocalFiniteElement,
+ MixedLocalGFEADOLCStiffness<DeformationFEBasis,
RealTuple<double,3>,
- OrientationFEBasis::LocalFiniteElement,
+ OrientationFEBasis,
Rotation<double,3> > localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
MixedGFEAssembler<DeformationFEBasis,
@@ -439,10 +447,10 @@ int main (int argc, char *argv[]) try
PythonFunction<FieldVector<double,dim>, FieldMatrix<double,3,3> > orientationDirichletValues(main.get("orientationDirichletValues"));
std::vector<FieldVector<double,3> > ddV;
- Functions::interpolate(deformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
+ ::Functions::interpolate(fufemDeformationFEBasis, ddV, deformationDirichletValues, deformationDirichletDofs);
std::vector<FieldMatrix<double,3,3> > dOV;
- Functions::interpolate(orientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
+ ::Functions::interpolate(fufemOrientationFEBasis, dOV, orientationDirichletValues, orientationDirichletDofs);
for (size_t j=0; j<xDisp.size(); j++)
if (deformationDirichletNodes[j][0])
@@ -464,8 +472,8 @@ int main (int argc, char *argv[]) try
// Output result of each homotopy step
std::stringstream iAsAscii;
iAsAscii << i+1;
- CosseratVTKWriter<GridType>::writeMixed<DeformationFEBasis,OrientationFEBasis>(deformationFEBasis,xDisp,
- orientationFEBasis,xOrient,
+ CosseratVTKWriter<GridType>::writeMixed<FufemDeformationFEBasis,FufemOrientationFEBasis>(fufemDeformationFEBasis,xDisp,
+ fufemOrientationFEBasis,xOrient,
resultPath + "mixed-cosserat_homotopy_" + iAsAscii.str());
}
diff --git a/src/rod-eoc.cc b/src/rod-eoc.cc
index 7cf35077..de8a21bd 100644
--- a/src/rod-eoc.cc
+++ b/src/rod-eoc.cc
@@ -1,5 +1,7 @@
#include <config.h>
+#define HAVE_CSTDDEF
+
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
@@ -8,6 +10,8 @@
#include <dune/istl/io.hh>
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
+
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
@@ -54,6 +58,14 @@ void solve (const GridType& grid,
const double E = parameters.get<double>("E");
const double nu = parameters.get<double>("nu");
+ // Create a function space basis
+ typedef Dune::Functions::PQKNodalBasis<typename GridType::LeafGridView, 1> FEBasis;
+ FEBasis feBasis(grid.leafGridView());
+
+ // Transitional: the same basis as a dune-fufem object
+ typedef DuneFunctionsBasis<FEBasis> FufemFEBasis;
+ FufemFEBasis fufemFeBasis(feBasis);
+
// Create a local assembler
RodLocalStiffness<OneDGrid::LeafGridView,double> localStiffness(grid.leafGridView(),
A, J1, J2, E, nu);
@@ -87,9 +99,9 @@ void solve (const GridType& grid,
// Create a solver for the rod problem
// ///////////////////////////////////////////
- RodAssembler<GridType::LeafGridView,3> rodAssembler(grid.leafGridView(), &localStiffness);
+ RodAssembler<FEBasis,3> rodAssembler(feBasis, &localStiffness);
- RiemannianTrustRegionSolver<GridType,RigidBodyMotion<double,3> > rodSolver;
+ RiemannianTrustRegionSolver<FEBasis,RigidBodyMotion<double,3> > rodSolver;
#if 1
rodSolver.setup(grid,
&rodAssembler,
diff --git a/src/rodobstacle.cc b/src/rodobstacle.cc
index 1cfa9893..a99625fc 100644
--- a/src/rodobstacle.cc
+++ b/src/rodobstacle.cc
@@ -1,5 +1,7 @@
#include <config.h>
+#define HAVE_CSTDDEF
+
#include <dune/common/bitsetvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
@@ -45,9 +47,9 @@ void setTrustRegionObstacles(double trustRegionRadius,
(trueObstacles[j].lower(k) < -1e10)
? std::min(-trustRegionRadius, trueObstacles[j].upper(k) - trustRegionRadius)
: trueObstacles[j].lower(k);
-
+
trustRegionObstacles[j].upper(k) =
- (trueObstacles[j].upper(k) > 1e10)
+ (trueObstacles[j].upper(k) > 1e10)
? std::max(trustRegionRadius,trueObstacles[j].lower(k) + trustRegionRadius)
: trueObstacles[j].upper(k);
@@ -88,10 +90,10 @@ int main (int argc, char *argv[]) try
const int baseIt = parameterSet.get<int>("baseIt");
const double tolerance = parameterSet.get<double>("tolerance");
const double baseTolerance = parameterSet.get<double>("baseTolerance");
-
+
// Problem settings
const int numRodBaseElements = parameterSet.get<int>("numRodBaseElements");
-
+
// ///////////////////////////////////////
// Create the two grids
// ///////////////////////////////////////
@@ -101,7 +103,7 @@ int main (int argc, char *argv[]) try
grid.globalRefine(minLevel);
std::vector<std::vector<BoxConstraint<double,3> > > trustRegionObstacles(minLevel+1);
- std::vector<BitSetVector<1> > hasObstacle(minLevel+1);
+ std::vector<BitSetVector<3> > hasObstacle(minLevel+1);
BitSetVector<blocksize> dirichletNodes;
// ////////////////////////////////
@@ -123,7 +125,7 @@ int main (int argc, char *argv[]) try
ProjectedBlockGSStep<MatrixType, CorrectionType> presmoother;
ProjectedBlockGSStep<MatrixType, CorrectionType> postsmoother;
- MonotoneMGStep<MatrixType, CorrectionType> multigridStep(1);
+ MonotoneMGStep<MatrixType, CorrectionType> multigridStep;
multigridStep.setMGType(mu, nu1, nu2);
multigridStep.ignoreNodes_ = &dirichletNodes;
@@ -163,16 +165,16 @@ int main (int argc, char *argv[]) try
// /////////////////////////////////////////////////////////////////////
// Refinement Loop
// /////////////////////////////////////////////////////////////////////
-
+
for (int toplevel=minLevel; toplevel<=maxLevel; toplevel++) {
-
+
std::cout << "####################################################" << std::endl;
std::cout << " Solving on level: " << toplevel << std::endl;
std::cout << "####################################################" << std::endl;
-
+
dirichletNodes.resize( grid.size(1) );
dirichletNodes.unsetAll();
-
+
dirichletNodes[0] = true;
dirichletNodes.back() = true;
@@ -183,41 +185,41 @@ int main (int argc, char *argv[]) try
MatrixType hessianMatrix;
RodAssembler<GridType::LeafGridView,2> rodAssembler(grid.leafGridView());
-
+
rodAssembler.setParameters(1, 350000, 350000);
-
+
MatrixIndexSet indices(grid.size(toplevel,1), grid.size(toplevel,1));
rodAssembler.getNeighborsPerVertex(indices);
indices.exportIdx(hessianMatrix);
-
+
rhs.resize(grid.size(toplevel,1));
corr.resize(grid.size(toplevel,1));
-
+
// //////////////////////////////////////////////////////////
// Create obstacles
// //////////////////////////////////////////////////////////
-
+
hasObstacle.resize(toplevel+1);
for (int i=0; i<hasObstacle.size(); i++) {
hasObstacle[i].resize(grid.size(i, 1));
hasObstacle[i].setAll();
}
-
+
std::vector<std::vector<BoxConstraint<double,3> > > trueObstacles(toplevel+1);
trustRegionObstacles.resize(toplevel+1);
-
+
for (int i=0; i<toplevel+1; i++) {
trueObstacles[i].resize(grid.size(i,1));
trustRegionObstacles[i].resize(grid.size(i,1));
}
-
+
for (int i=0; i<trueObstacles[toplevel].size(); i++) {
trueObstacles[toplevel][i].clear();
//trueObstacles[toplevel][i].val[0] = - x[i][0];
trueObstacles[toplevel][i].upper(0) = 0.1 - x[i].r[0];
}
-
+
trustRegionObstacles.resize(toplevel+1);
for (int i=0; i<=toplevel; i++)
@@ -248,17 +250,17 @@ int main (int argc, char *argv[]) try
for (int i=0; i<maxNewtonSteps; i++) {
std::cout << "-----------------------------------------------------------------------------" << std::endl;
- std::cout << " Trust-Region Step Number: " << i
+ std::cout << " Trust-Region Step Number: " << i
<< ", radius: " << trustRegionRadius
<< ", energy: " << rodAssembler.computeEnergy(x) << std::endl;
std::cout << "-----------------------------------------------------------------------------" << std::endl;
rhs = 0;
corr = 0;
-
+
rodAssembler.assembleGradient(x, rhs);
rodAssembler.assembleMatrix(x, hessianMatrix);
-
+
rhs *= -1;
// Create trust-region obstacle on grid0.maxLevel()
@@ -292,18 +294,18 @@ int main (int argc, char *argv[]) try
// ////////////////////////////////////////////////////
SolutionType newIterate = x;
- for (int j=0; j<newIterate.size(); j++)
+ for (int j=0; j<newIterate.size(); j++)
newIterate[j] = RigidBodyMotion<double,2>::exp(newIterate[j], corr[j]);
/** \todo Don't always recompute oldEnergy */
- double oldEnergy = rodAssembler.computeEnergy(x);
- double energy = rodAssembler.computeEnergy(newIterate);
+ double oldEnergy = rodAssembler.computeEnergy(x);
+ double energy = rodAssembler.computeEnergy(newIterate);
- if (energy >= oldEnergy)
+ if (energy >= oldEnergy)
DUNE_THROW(SolverError, "Richtung ist keine Abstiegsrichtung!");
-
+
// Add correction to the current solution
- for (int j=0; j<x.size(); j++)
+ for (int j=0; j<x.size(); j++)
x[j] = RigidBodyMotion<double,2>::exp(x[j], corr[j]);
// Subtract correction from the current obstacle
@@ -311,31 +313,31 @@ int main (int argc, char *argv[]) try
trueObstacles[grid.maxLevel()][k] -= corr[k];
}
-
+
// //////////////////////////////
// Output result
// //////////////////////////////
-
+
// Write Lagrange multiplyers
std::stringstream levelAsAscii;
levelAsAscii << toplevel;
std::string lagrangeFilename = "pressure/lagrange_" + levelAsAscii.str();
std::ofstream lagrangeFile(lagrangeFilename.c_str());
-
+
CorrectionType lagrangeMultipliers;
rodAssembler.assembleGradient(x, lagrangeMultipliers);
lagrangeFile << lagrangeMultipliers << std::endl;
-
+
// Write result grid
std::string solutionFilename = "solutions/rod_" + levelAsAscii.str() + ".result";
writeRod(x, solutionFilename);
-
+
// ////////////////////////////////////////////////////////////////////////////
// Refine locally and transfer the current solution to the new leaf level
// ////////////////////////////////////////////////////////////////////////////
-
+
GeometricEstimator<GridType> estimator;
-
+
estimator.estimate(grid, (toplevel<=minLevel) ? refineAll : refineCondition);
std::cout << " #### WARNING: function not transferred to the next level! #### " << std::endl;
--
GitLab