From e57a9962fb4e88a2d9d5bfb860b7dc6fcff346a7 Mon Sep 17 00:00:00 2001 From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de> Date: Wed, 25 Mar 2020 14:05:37 +0100 Subject: [PATCH] Add option for adaptive refinement, write the grid using an unstructured grid writer and create a file containing the final deformation The grid and the file can then be read by dune-gfe Make dune-vtk optional --- dune.module | 2 +- .../finite-strain-elasticity-bending.parset | 53 ++++++++----- src/finite-strain-elasticity.cc | 78 +++++++++++++++---- 3 files changed, 97 insertions(+), 36 deletions(-) diff --git a/dune.module b/dune.module index 1ba8812..1029d09 100644 --- a/dune.module +++ b/dune.module @@ -8,4 +8,4 @@ Version: 2.7-git Maintainer: oliver.sander@tu-dresden.de, youett@mi.fu-berlin.de #depending on Depends:dune-common dune-grid dune-istl dune-solvers dune-fufem dune-geometry dune-functions -Suggests: dune-uggrid dune-parmg +Suggests: dune-uggrid dune-parmg dune-vtk diff --git a/problems/finite-strain-elasticity-bending.parset b/problems/finite-strain-elasticity-bending.parset index 8b0d639..5e9b930 100644 --- a/problems/finite-strain-elasticity-bending.parset +++ b/problems/finite-strain-elasticity-bending.parset @@ -3,14 +3,21 @@ ############################################# structuredGrid = true + +# bounding box lower = 0 0 0 -upper = 5 1 1 -elements = 20 5 5 +upper = 200 100 200 + +elements = 10 5 5 + +# Number of grid levels: refinement levels of the surface shell part, refined using hierarchic refinement +numLevels = 2 + #Overlap indicator for the grid when doing parallel calculations; when set to false, all communication will be done assuming there is no overlap when assembling the matrix and the rhs overlap = false -# Number of grid levels -numLevels = 1 +# Adaptive Refinement +adaptiveRefinement = true ############################################# # Solver parameters @@ -23,13 +30,13 @@ numHomotopySteps = 1 tolerance = 1e-6 # Max number of steps of the trust region solver -maxTrustRegionSteps = 500 +maxTrustRegionSteps = 200 # Initial trust-region radius -initialTrustRegionRadius = 0.1 +initialTrustRegionRadius = 10 # Number of multigrid iterations per trust-region step -numIt = 1000 +numIt = 100 # Number of presmoothing steps nu1 = 3 @@ -56,31 +63,31 @@ baseTolerance = 1e-8 # Material parameters ############################ -energy = stvenantkirchhoff +energy = mooneyrivlin ## For the Wriggers L-shape example [materialParameters] -# Lame parameters +# Lame parameters for stvenantkirchhoff # corresponds to E = 71240 N/mm^2, nu=0.31 # However, we use units N/m^2 -mu = 2.7191e+6 -lambda = 4.4364e+6 +mu = 2.7191e+4 +lambda = 4.4364e+4 mooneyrivlin_a = 2.7191e+6 mooneyrivlin_b = 2.7191e+6 mooneyrivlin_c = 2.7191e+6 -mooneyrivlin_10 = -7.28e+5 -mooneyrivlin_01 = 9.17e+5 -mooneyrivlin_20 = 1.23e+5 -mooneyrivlin_02 = 8.15e+5 -mooneyrivlin_11 = -5.14e+5 +mooneyrivlin_10 = -7.28e+4 +mooneyrivlin_01 = 9.17e+4 +mooneyrivlin_20 = 1.23e+4 +mooneyrivlin_02 = 8.15e+4 +mooneyrivlin_11 = -5.14e+4 mooneyrivlin_30 = 0 mooneyrivlin_21 = 0 mooneyrivlin_12 = 0 mooneyrivlin_03 = 0 -mooneyrivlin_k = 1e+6 # volume-preserving parameter +mooneyrivlin_k = 1e+5 # volume-preserving parameter mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute the Mooney-Rivlin-Energy [] @@ -95,12 +102,16 @@ problem = wriggers-l-shape # x is the vertex coordinate dirichletVerticesPredicate = "x[0] < 0.01" -### Python predicate specifying all Dirichlet grid vertices +### Python predicate specifying all surface shell grid vertices +# x is the vertex coordinate +adaptiveRefinementPredicate = "x[2] > 199.99" + +### Python predicate specifying all Neumann grid vertices # x is the vertex coordinate -neumannVerticesPredicate = "x[0] > 4.99" +neumannVerticesPredicate = "x[0] > 199.99" ### Neumann values, if needed -neumannValues = 0 5e4 0 +neumannValues = 15e3 0 0 # Initial deformation -initialDeformation = "[x[0], x[1], x[2]]" +initialDeformation = "x" diff --git a/src/finite-strain-elasticity.cc b/src/finite-strain-elasticity.cc index ae57f84..2cf3c74 100644 --- a/src/finite-strain-elasticity.cc +++ b/src/finite-strain-elasticity.cc @@ -42,6 +42,13 @@ #include <dune/elasticity/materials/neumannenergy.hh> #include <dune/elasticity/materials/sumenergy.hh> +#if HAVE_DUNE_VTK +#include <dune/vtk/writers/vtkunstructuredgridwriter.hh> +#endif + +#include <iostream> +#include <fstream> + // grid dimension const int dim = 3; @@ -97,7 +104,7 @@ int main (int argc, char *argv[]) try ParameterTreeParser::readOptions(argc, argv, parameterSet); // read solver settings - const int numLevels = parameterSet.get<int>("numLevels"); + 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"); @@ -121,6 +128,8 @@ int main (int argc, char *argv[]) try std::shared_ptr<GridType> grid; + const bool adaptiveRefinement = parameterSet.get<bool>("adaptiveRefinement"); + FieldVector<double,dim> lower(0), upper(1); if (parameterSet.get<bool>("structuredGrid")) { @@ -129,17 +138,49 @@ int main (int argc, char *argv[]) try upper = parameterSet.get<FieldVector<double,dim> >("upper"); std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements"); - grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, 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 = std::shared_ptr<GridType>(GmshReader<GridType>::read(path + "/" + gridFile)); } - grid->globalRefine(numLevels-1); + // 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)); + + // 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)); + + if (adaptiveRefinement) { + // Same for the boundary that will get adaptively refined + lambda = std::string("lambda x: (") + parameterSet.get<std::string>("adaptiveRefinementPredicate", "0") + std::string(")"); + PythonFunction<FieldVector<double,dim>, bool> pythonSurfaceShellVertices(Python::evaluate(lambda)); + + grid->setRefinementType(GridType::RefinementType::COPY); + + while (numLevels > 0) { + for (auto&& e : elements(grid->leafGridView())){ + bool isSurfaceShell = false; + for (int i = 0; i < e.geometry().corners(); i++) { + isSurfaceShell = isSurfaceShell || pythonSurfaceShellVertices(e.geometry().corner(i)); + } + grid->mark(isSurfaceShell ? 1 : 0,e); + } + + grid->adapt(); - grid->loadBalance(); + grid->loadBalance(); + + numLevels--; + } + } else { + grid->globalRefine(numLevels-1); + grid->loadBalance(); + } if (mpiHelper.rank()==0) std::cout << "There are " << grid->leafGridView().comm().size() << " processes" << std::endl; @@ -165,15 +206,6 @@ int main (int argc, char *argv[]) try const GridView::IndexSet& indexSet = gridView.indexSet(); - // 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)); - - // 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)); - for (auto&& v : vertices(gridView)) { bool isDirichlet; @@ -384,7 +416,25 @@ int main (int argc, char *argv[]) try vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim)); vtkWriter.write(resultPath + "finite-strain_homotopy_" + std::to_string(i+1)); } - +#if HAVE_DUNE_VTK + if (grid->leafGridView().comm().size() == 1) { + displacement = x; + displacement -= identity; + auto finalDisplacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement); + //Write out the grid using an unstructured GridWriter + VtkUnstructuredGridWriter<GridView> vtkUnstructuredGridWriter(gridView, Vtk::ASCII); + vtkUnstructuredGridWriter.addPointData(finalDisplacementFunction, "displacement"); + vtkUnstructuredGridWriter.write(resultPath + "finite-strain_homotopy_original"); + std::ofstream file; + file.open("deformation"); + for (auto& d : displacement){ + file << d << "\n"; + } + file.close(); + } else if (mpiHelper.rank() == 0) { + std::cout << "Not writing out the deformation to a file, this currently only works for sequential calculations" << std::endl; + } +#endif if (mpiHelper.rank()==0) std::cout << "Complete duration: " << homotopyTimer.elapsed() << " sec." << std::endl; -- GitLab