diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc
index 681fbc7b73c3e2882d96d17fe1512385813f6c65..b6471e4937716cd46bd4ad3bd9f0137d455e1499 100644
--- a/src/finite-strain-elasticity.cc
+++ b/src/finite-strain-elasticity.cc
@@ -46,194 +46,194 @@ using namespace Dune;
struct NeumannFunction
: public Dune::VirtualFunction<FieldVector<double,dim>, FieldVector<double,3> >
{
- NeumannFunction(const FieldVector<double,3> values,
- double homotopyParameter)
- : values_(values),
- homotopyParameter_(homotopyParameter)
- {}
-
- void evaluate(const FieldVector<double, dim>& x, FieldVector<double,3>& out) const {
- out = 0;
- out.axpy(-homotopyParameter_, values_);
- }
-
- FieldVector<double,3> values_;
- double homotopyParameter_;
+ NeumannFunction(const FieldVector<double,3> values,
+ double homotopyParameter)
+ : values_(values),
+ homotopyParameter_(homotopyParameter)
+ {}
+
+ void evaluate(const FieldVector<double, dim>& x, FieldVector<double,3>& out) const
+ {
+ out = 0;
+ out.axpy(-homotopyParameter_, values_);
+ }
+
+ FieldVector<double,3> values_;
+ double homotopyParameter_;
};
int main (int argc, char *argv[]) try
{
- // initialize MPI, finalize is done automatically on exit
- Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+ // initialize MPI, finalize is done automatically on exit
+ Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
- // Start Python interpreter
- Python::start();
- Python::Reference main = Python::import("__main__");
- Python::run("import math");
+ // Start Python interpreter
+ Python::start();
+ Python::Reference main = Python::import("__main__");
+ Python::run("import math");
- //feenableexcept(FE_INVALID);
- Python::runStream()
+ //feenableexcept(FE_INVALID);
+ Python::runStream()
<< std::endl << "import sys"
<< std::endl << "sys.path.append('/home/sander/dune/dune-gfe/')"
<< std::endl;
- typedef std::vector<FieldVector<double,dim> > SolutionType;
+ typedef std::vector<FieldVector<double,dim> > SolutionType;
- // parse data file
- ParameterTree parameterSet;
+ // parse data file
+ ParameterTree parameterSet;
// if (argc != 2)
// DUNE_THROW(Exception, "Usage: ./hencky-material <parameter file>");
- ParameterTreeParser::readINITree(argv[1], parameterSet);
+ ParameterTreeParser::readINITree(argv[1], parameterSet);
- ParameterTreeParser::readOptions(argc, argv, parameterSet);
+ ParameterTreeParser::readOptions(argc, argv, parameterSet);
- // read solver settings
- const int numLevels = parameterSet.get<int>("numLevels");
- int numHomotopySteps = parameterSet.get<int>("numHomotopySteps");
- const double tolerance = parameterSet.get<double>("tolerance");
- const int maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
- const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
- const int multigridIterations = parameterSet.get<int>("numIt");
- 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");
- std::string resultPath = parameterSet.get("resultPath", "");
+ // read solver settings
+ const int numLevels = parameterSet.get<int>("numLevels");
+ int numHomotopySteps = parameterSet.get<int>("numHomotopySteps");
+ const double tolerance = parameterSet.get<double>("tolerance");
+ const int maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
+ const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+ const int multigridIterations = parameterSet.get<int>("numIt");
+ 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");
+ std::string resultPath = parameterSet.get("resultPath", "");
- // ///////////////////////////////////////
- // Create the grid
- // ///////////////////////////////////////
- typedef UGGrid<dim> GridType;
+ // ///////////////////////////////////////
+ // Create the grid
+ // ///////////////////////////////////////
+ typedef UGGrid<dim> GridType;
- shared_ptr<GridType> grid;
+ shared_ptr<GridType> grid;
- FieldVector<double,dim> lower(0), upper(1);
+ FieldVector<double,dim> lower(0), upper(1);
- if (parameterSet.get<bool>("structuredGrid")) {
+ if (parameterSet.get<bool>("structuredGrid")) {
- lower = parameterSet.get<FieldVector<double,dim> >("lower");
- upper = parameterSet.get<FieldVector<double,dim> >("upper");
+ lower = parameterSet.get<FieldVector<double,dim> >("lower");
+ upper = parameterSet.get<FieldVector<double,dim> >("upper");
- array<unsigned int,dim> elements = parameterSet.get<array<unsigned int,dim> >("elements");
- grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+ array<unsigned int,dim> elements = parameterSet.get<array<unsigned int,dim> >("elements");
+ grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
- } else {
- std::string path = parameterSet.get<std::string>("path");
- std::string gridFile = parameterSet.get<std::string>("gridFile");
- grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
- }
+ } else {
+ std::string path = parameterSet.get<std::string>("path");
+ std::string gridFile = parameterSet.get<std::string>("gridFile");
+ grid = shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile));
+ }
- grid->globalRefine(numLevels-1);
+ grid->globalRefine(numLevels-1);
- grid->loadBalance();
+ grid->loadBalance();
- if (mpiHelper.rank()==0)
- std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
+ if (mpiHelper.rank()==0)
+ std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl;
- typedef GridType::LeafGridView GridView;
- GridView gridView = grid->leafGridView();
+ typedef GridType::LeafGridView GridView;
+ GridView gridView = grid->leafGridView();
// typedef P1NodalBasis<GridView,double> FEBasis;
- typedef P1NodalBasis<GridView,double> FEBasis;
- FEBasis feBasis(gridView);
-
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
+ typedef P1NodalBasis<GridView,double> FEBasis;
+ FEBasis feBasis(gridView);
- BitSetVector<1> dirichletVertices(gridView.size(dim), false);
- BitSetVector<1> neumannVertices(gridView.size(dim), false);
+ // /////////////////////////////////////////
+ // Read Dirichlet values
+ // /////////////////////////////////////////
- GridType::Codim<dim>::LeafIterator vIt = gridView.begin<dim>();
- GridType::Codim<dim>::LeafIterator vEndIt = gridView.end<dim>();
+ BitSetVector<1> dirichletVertices(gridView.size(dim), false);
+ BitSetVector<1> neumannVertices(gridView.size(dim), false);
- const GridView::IndexSet& indexSet = gridView.indexSet();
+ GridType::Codim<dim>::LeafIterator vIt = gridView.begin<dim>();
+ GridType::Codim<dim>::LeafIterator vEndIt = gridView.end<dim>();
- // Make Python function that computes which vertices are on the Dirichlet boundary,
- // based on the vertex positions.
- std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
- PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+ const GridView::IndexSet& indexSet = gridView.indexSet();
- // Same for the Neumann boundary
- lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
- PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
+ // Make Python function that computes which vertices are on the Dirichlet boundary,
+ // based on the vertex positions.
+ std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
- for (; vIt!=vEndIt; ++vIt) {
+ // Same for the Neumann boundary
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
- bool isDirichlet;
- pythonDirichletVertices.evaluate(vIt->geometry().corner(0), isDirichlet);
- dirichletVertices[indexSet.index(*vIt)] = isDirichlet;
+ for (; vIt!=vEndIt; ++vIt)
+ {
+ bool isDirichlet;
+ pythonDirichletVertices.evaluate(vIt->geometry().corner(0), isDirichlet);
+ dirichletVertices[indexSet.index(*vIt)] = isDirichlet;
- bool isNeumann;
- pythonNeumannVertices.evaluate(vIt->geometry().corner(0), isNeumann);
- neumannVertices[indexSet.index(*vIt)] = isNeumann;
+ bool isNeumann;
+ pythonNeumannVertices.evaluate(vIt->geometry().corner(0), isNeumann);
+ neumannVertices[indexSet.index(*vIt)] = isNeumann;
+ }
- }
+ BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
+ BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
- BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
- BoundaryPatch<GridView> neumannBoundary(gridView, neumannVertices);
+ if (mpiHelper.rank()==0)
+ std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
- if (mpiHelper.rank()==0)
- std::cout << "Neumann boundary has " << neumannBoundary.numFaces() << " faces\n";
+ BitSetVector<1> dirichletNodes(feBasis.size(), false);
+ constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
- BitSetVector<1> dirichletNodes(feBasis.size(), false);
- constructBoundaryDofs(dirichletBoundary,feBasis,dirichletNodes);
+ BitSetVector<1> neumannNodes(feBasis.size(), false);
+ constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
- BitSetVector<1> neumannNodes(feBasis.size(), false);
- constructBoundaryDofs(neumannBoundary,feBasis,neumannNodes);
+ BitSetVector<dim> dirichletDofs(feBasis.size(), false);
+ for (size_t i=0; i<feBasis.size(); i++)
+ if (dirichletNodes[i][0])
+ for (int j=0; j<dim; j++)
+ dirichletDofs[i][j] = true;
- BitSetVector<dim> dirichletDofs(feBasis.size(), false);
- for (size_t i=0; i<feBasis.size(); i++)
- if (dirichletNodes[i][0])
- for (int j=0; j<dim; j++)
- dirichletDofs[i][j] = true;
+ // //////////////////////////
+ // Initial iterate
+ // //////////////////////////
- // //////////////////////////
- // Initial iterate
- // //////////////////////////
+ SolutionType x(feBasis.size());
- SolutionType x(feBasis.size());
+ lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > pythonInitialDeformation(Python::evaluate(lambda));
- 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(feBasis, v, pythonInitialDeformation);
- std::vector<FieldVector<double,3> > v;
- Functions::interpolate(feBasis, v, pythonInitialDeformation);
+ for (size_t i=0; i<x.size(); i++)
+ x[i] = v[i];
- for (size_t i=0; i<x.size(); i++)
- x[i] = v[i];
+ ////////////////////////////////////////////////////////
+ // Main homotopy loop
+ ////////////////////////////////////////////////////////
- ////////////////////////////////////////////////////////
- // Main homotopy loop
- ////////////////////////////////////////////////////////
+ // Output initial iterate (of homotopy loop)
+ VTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView());
+ BlockVector<FieldVector<double,3> > displacement(x.size());
+ for (auto it = grid->leafGridView().template begin<dim>(); it != grid->leafGridView().template end<dim>(); ++it)
+ {
+ size_t idx = grid->leafGridView().indexSet().index(*it);
+ displacement[idx] = x[idx] - it->geometry().corner(0);
- // Output initial iterate (of homotopy loop)
- VTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView());
- BlockVector<FieldVector<double,3> > displacement(x.size());
- for (auto it = grid->leafGridView().template begin<dim>(); it != grid->leafGridView().template end<dim>(); ++it)
- {
- size_t idx = grid->leafGridView().indexSet().index(*it);
- displacement[idx] = x[idx] - it->geometry().corner(0);
-
- //std::cout << "idx: " << idx << " coordinate: " << it->geometry().corner(0) << std::endl;
- }
-
- Dune::shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > > vtkDisplacement
- = Dune::make_shared<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > >
- (feBasis, displacement, "Displacement");
- vtkWriter.addVertexData(vtkDisplacement);
- vtkWriter.write(resultPath + "hencky_homotopy_0");
+ //std::cout << "idx: " << idx << " coordinate: " << it->geometry().corner(0) << std::endl;
+ }
- for (int i=0; i<numHomotopySteps; i++) {
+ Dune::shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > > vtkDisplacement
+ = Dune::make_shared<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > >
+ (feBasis, displacement, "Displacement");
+ vtkWriter.addVertexData(vtkDisplacement);
+ vtkWriter.write(resultPath + "hencky_homotopy_0");
- double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
- if (mpiHelper.rank()==0)
- std::cout << "Homotopy step: " << i << ", parameter: " << homotopyParameter << std::endl;
+ for (int i=0; i<numHomotopySteps; i++)
+ {
+ double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
+ if (mpiHelper.rank()==0)
+ std::cout << "Homotopy step: " << i << ", parameter: " << homotopyParameter << std::endl;
// ////////////////////////////////////////////////////////////
@@ -243,23 +243,24 @@ int main (int argc, char *argv[]) try
const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
shared_ptr<NeumannFunction> neumannFunction;
if (parameterSet.hasKey("neumannValues"))
- neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
- homotopyParameter);
+ neumannFunction = make_shared<NeumannFunction>(parameterSet.get<FieldVector<double,3> >("neumannValues"),
+ homotopyParameter);
- std::cout << "Neumann values: " << parameterSet.get<FieldVector<double,3> >("neumannValues") << std::endl;
+ std::cout << "Neumann values: " << parameterSet.get<FieldVector<double,3> >("neumannValues") << std::endl;
- if (mpiHelper.rank() == 0) {
- std::cout << "Material parameters:" << std::endl;
- materialParameters.report();
- }
+ if (mpiHelper.rank() == 0)
+ {
+ std::cout << "Material parameters:" << std::endl;
+ materialParameters.report();
+ }
// Assembler using ADOL-C
StVenantKirchhoffEnergy<GridView,
- FEBasis::LocalFiniteElement,
- adouble> henckyEnergy(materialParameters,
- &neumannBoundary,
- neumannFunction.get());
+ FEBasis::LocalFiniteElement,
+ adouble> henckyEnergy(materialParameters,
+ &neumannBoundary,
+ neumannFunction.get());
LocalADOLCStiffness<GridView,
FEBasis::LocalFiniteElement,
@@ -286,59 +287,57 @@ int main (int argc, char *argv[]) try
baseTolerance
);
- ////////////////////////////////////////////////////////
- // Set Dirichlet values
- ////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////
+ // Set Dirichlet values
+ ////////////////////////////////////////////////////////
- Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
+ Python::Reference dirichletValuesClass = Python::import(parameterSet.get<std::string>("problem") + "-dirichlet-values");
- Python::Callable C = dirichletValuesClass.get("DirichletValues");
+ Python::Callable C = dirichletValuesClass.get("DirichletValues");
- // Call a constructor.
- Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
+ // Call a constructor.
+ Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
- // Extract object member functions as Dune functions
- PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > dirichletValues(dirichletValuesPythonObject.get("dirichletValues"));
+ // Extract object member functions as Dune functions
+ PythonFunction<FieldVector<double,dim>, FieldVector<double,3> > dirichletValues(dirichletValuesPythonObject.get("dirichletValues"));
- std::vector<FieldVector<double,3> > ddV;
- Functions::interpolate(feBasis, ddV, dirichletValues, dirichletDofs);
+ std::vector<FieldVector<double,3> > ddV;
+ Functions::interpolate(feBasis, ddV, dirichletValues, dirichletDofs);
- for (size_t j=0; j<x.size(); j++)
- if (dirichletNodes[j][0])
- x[j] = ddV[j];
+ for (size_t j=0; j<x.size(); j++)
+ if (dirichletNodes[j][0])
+ x[j] = ddV[j];
- // /////////////////////////////////////////////////////
- // Solve!
- // /////////////////////////////////////////////////////
+ // /////////////////////////////////////////////////////
+ // Solve!
+ // /////////////////////////////////////////////////////
- solver.setInitialIterate(x);
- solver.solve();
+ solver.setInitialIterate(x);
+ solver.solve();
- x = solver.getSol();
+ x = solver.getSol();
- // Output result of each homotopy step
- VTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView());
- BlockVector<FieldVector<double,3> > displacement(x.size());
- for (auto it = grid->leafGridView().template begin<dim>(); it != grid->leafGridView().template end<dim>(); ++it)
- {
- size_t idx = grid->leafGridView().indexSet().index(*it);
- displacement[idx] = x[idx] - it->geometry().corner(0);
- }
+ // Output result of each homotopy step
+ VTKWriter<GridType::LeafGridView> vtkWriter(grid->leafGridView());
+ BlockVector<FieldVector<double,3> > displacement(x.size());
+ for (auto it = grid->leafGridView().template begin<dim>(); it != grid->leafGridView().template end<dim>(); ++it)
+ {
+ size_t idx = grid->leafGridView().indexSet().index(*it);
+ displacement[idx] = x[idx] - it->geometry().corner(0);
+ }
- /////////////////////////////////
- // Output result
- /////////////////////////////////
+ /////////////////////////////////
+ // Output result
+ /////////////////////////////////
- Dune::shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > > vtkDisplacement
- = Dune::make_shared<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > >
+ Dune::shared_ptr<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > > vtkDisplacement
+ = Dune::make_shared<VTKBasisGridFunction<FEBasis,BlockVector<FieldVector<double,3> > > >
(feBasis, displacement, "Displacement");
- vtkWriter.addVertexData(vtkDisplacement);
- vtkWriter.write(resultPath + "hencky_homotopy_" + std::to_string(i+1));
-
- }
+ vtkWriter.addVertexData(vtkDisplacement);
+ vtkWriter.write(resultPath + "hencky_homotopy_" + std::to_string(i+1));
- } catch (Exception e) {
+ }
+} catch (Exception e) {
std::cout << e << std::endl;
-
- }
+}