Skip to content
Snippets Groups Projects
Commit e90d6953 authored by lisa_julia.nebel_at_tu-dresden.de's avatar lisa_julia.nebel_at_tu-dresden.de
Browse files

Add option for adaptive refinement, write the grid using an unstructured grid...

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
parent ca987c8d
No related branches found
No related tags found
No related merge requests found
......@@ -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
......@@ -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"
......@@ -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,18 +138,50 @@ 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();
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;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment