Commit 364f6038 authored by Peter Bastian's avatar Peter Bastian
Browse files

integrated geneo setup into new example

parent 5fcf41a9
################################
# Dune module information file #
################################
......@@ -15,6 +7,6 @@ Module: dune-fiber-elasticity
Version: 2.7
Maintainer: Dominic Kempf
# Required build dependencies
Depends: dune-pdelab
Depends: dune-pdelab dune-ftworth
# Optional build dependencies
......@@ -70,7 +70,7 @@ class FibreReinforcedBulkOperator
virtual void alpha_skeleton(const IG& ig, const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s, const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n, R& r_s, R& r_n) const
{
bulkoperator.alpha_skeleton(ig, lfsu_s, x_s, lfsv_s, lfsu_n, x_n, lfsv_n, r_s, r_n);
fibreoperator.alpha_skeleton(ig, lfsu_s, x_s, lfsv_s, lfsu_n, x_n, lfsv_n, r_s, r_n);
//fibreoperator.alpha_skeleton(ig, lfsu_s, x_s, lfsv_s, lfsu_n, x_n, lfsv_n, r_s, r_n);
}
private:
......@@ -78,4 +78,4 @@ class FibreReinforcedBulkOperator
EulerBernoulli2DLocalOperator<GFS, FGFS> fibreoperator;
};
#endif
\ No newline at end of file
#endif
add_executable(example example.cc)
target_link_dune_default_libraries(example)
add_executable(example_geneo example_geneo.cc)
target_link_dune_default_libraries(example_geneo)
dune_symlink_to_source_files(FILES example.yml)
......@@ -6,8 +6,8 @@ grid:
- 4.0
- 1.0
N:
- 40
- 10
- 20
- 20
refinement: 0
material:
-
......@@ -36,3 +36,10 @@ reinforced_operator:
- 0.749
youngs_modulus: 10000
stabilization_parameter: 1000
geneo:
subdomains:
- 4
- 1
overlap: 2
extensionmethod: vertex
drop_small_overlap: 1
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
// dune-fiber-elasticity specific includes you probably want to add
#include<yaml-cpp/yaml.h>
#include<dune/fiber-elasticity/material.hh>
#include<dune/fiber-elasticity/operator.hh>
// Stuff that I use for setting up the problem
#include<dune/fiber-elasticity/callabletree.hh>
#include<dune/fiber-elasticity/vonmises.hh>
// General Dune stuff
#include<dune/common/parallel/mpihelper.hh>
#include<dune/grid/io/file/gmshreader.hh>
#include<dune/grid/uggrid.hh>
#include<dune/istl/preconditioners.hh>
#include<dune/pdelab.hh>
#include<memory>
#include<vector>
// ftworth includes
#include<dune/ftworth/partitioner.hh>
#include<dune/ftworth/assemblewrapper.hh>
#include<dune/ftworth/subdomainutilities.hh>
#include<dune/ftworth/cdediscretizer.hh>
#include<dune/ftworth/multilevel_geneo_preconditioner.hh>
int main(int argc, char** argv)
{
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
// Parse configuration
if (argc==1)
{
std::cout << "./example <yml file>" << std::endl;
exit(1);
}
YAML::Node config = YAML::LoadFile(argv[1]);
// Construct a 2D UG grid - only sequential here
using Grid = Dune::UGGrid<2>;
auto ll = config["grid"]["lowerleft"].template as<Dune::FieldVector<double, 2>>(Dune::FieldVector<double, 2>(0.0));
auto ur = config["grid"]["upperright"].template as<Dune::FieldVector<double, 2>>(Dune::FieldVector<double, 2>(1.0));
auto N = config["grid"]["N"].template as<std::array<unsigned int, 2>>(Dune::filledArray<2, unsigned int>(10));
auto grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(ll, ur, N);
grid->loadBalance();
grid->globalRefine(config["grid"]["refinement"].template as<int>(0));
auto physical = std::make_shared<std::vector<int>>(grid->size(0), 0);
auto gv = grid->leafGridView();
using EntitySet = Dune::PDELab::OverlappingEntitySet<Grid::Traits::LeafGridView>;
EntitySet es(gv);
// PDELab stuff
using FiniteElementMap = Dune::PDELab::PkLocalFiniteElementMap<EntitySet, double, double, 2>;
using ConstraintsAssembler = Dune::PDELab::OverlappingConformingDirichletConstraints;
using VectorBackend = Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::none>;
using LeafVectorBackend = Dune::PDELab::ISTL::VectorBackend<>;
using VectorGridFunctionSpace = Dune::PDELab::VectorGridFunctionSpace<EntitySet, FiniteElementMap, 2, VectorBackend, LeafVectorBackend, ConstraintsAssembler>;
using ConstraintsContainer = typename VectorGridFunctionSpace::template ConstraintsContainer<double>::Type;
FiniteElementMap fem(es);
VectorGridFunctionSpace vgfs(es, fem);
ConstraintsContainer cc;
// Parsing the material description
auto material = parse_material<double>(es, physical, config["material"], config);
// Interpolation of initial conditions
using Vector = Dune::PDELab::Backend::Vector<VectorGridFunctionSpace, double>;
Vector displacement(vgfs, 0.0);
Vector body_force(vgfs, 0.0);
Vector traction_force(vgfs, 0.0);
auto body_force_func = Dune::BlockLab::makeGridFunctionTreeFromCallables(vgfs,
[](auto, auto){ return 0; }, [](auto, auto){ return -1; });
Dune::PDELab::interpolate(body_force_func, vgfs, body_force);
auto cc_func = Dune::BlockLab::makeBoundaryConditionTreeFromCallables(vgfs,
[](auto e, auto x){ return e.geometry().global(x)[0] < 1e-8; },
[](auto e, auto x){ return e.geometry().global(x)[0] < 1e-8; });
Dune::PDELab::constraints(cc_func, vgfs, cc);
std::cout << "constrained dofs=" << cc.size() << " of " << vgfs.globalSize() << std::endl;
// Construction of operators
using LocalOperator = FibreReinforcedBulkOperator<VectorGridFunctionSpace, VectorGridFunctionSpace, VectorGridFunctionSpace, 2>;
LocalOperator lop(Dune::stackobject_to_shared_ptr(vgfs), config["reinforced_operator"], material);
Dune::PDELab::ISTL::BCRSMatrixBackend<> mb(36);
using GridOperator = Dune::PDELab::GridOperator<VectorGridFunctionSpace,
VectorGridFunctionSpace,
LocalOperator,
Dune::PDELab::ISTL::BCRSMatrixBackend<>,
double,
double,
double,
ConstraintsContainer,
ConstraintsContainer>;
GridOperator go(vgfs, cc, vgfs, cc, lop, mb);
lop.setCoefficientTraction(Dune::stackobject_to_shared_ptr(vgfs), Dune::stackobject_to_shared_ptr(traction_force));
lop.setCoefficientForce(Dune::stackobject_to_shared_ptr(vgfs), Dune::stackobject_to_shared_ptr(body_force));
// Linear Solver Setup
using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_UMFPack;
using StationaryLinearProblemSolver = Dune::PDELab::StationaryLinearProblemSolver<GridOperator, LinearSolver, Vector>;
LinearSolver ls(0);
StationaryLinearProblemSolver slp(go, ls, displacement, 1e-12);
//==================
// Begin Geneo stuff
//==================
// read parameters
std::vector<size_t> subdomains = config["geneo"]["subdomains"].template as<std::vector<size_t>>(std::vector<size_t>(2,1));
std::cout << "subdomains: "; for (auto p : subdomains) std::cout << p << " "; std::cout << std::endl;
size_t overlap = config["geneo"]["overlap"].template as<size_t>(1);
std::string extensionmethod = config["geneo"]["extensionmethod"].template as<std::string>(std::string("vertex"));
size_t drop_small_overlap = config["geneo"]["drop_small_overlap"].template as<size_t>(1);
// load balancing and subdomain information for geneo
auto pdd = std::make_shared<DomainDecomposition>(gv,Dune::MPIHelper::getLocalCommunicator(),subdomains[0],overlap,extensionmethod,drop_small_overlap,false);
auto pdddm = std::make_shared<DomainDecompositionDOFMapper>(vgfs,pdd);
// set up linear system in PDELab components
using Matrix = typename GridOperator::Jacobian;
using ISTLV = Dune::PDELab::Backend::Native<Vector>;
using ISTLM = Dune::PDELab::Backend::Native<Matrix>;
ISTLV& istldisplacement = Dune::PDELab::Backend::native(displacement);
Matrix A(go); // this makes the sparsity pattern
A = 0.0; // clear matrix
ISTLM& istlA = Dune::PDELab::Backend::native(A);
std::shared_ptr<ISTLM> pistlA = stackobject_to_shared_ptr(Dune::PDELab::Backend::native(istlA));
// residual
Vector residual(vgfs);
residual = 0.0;
go.residual(displacement,residual);
ISTLV& istlr = Dune::PDELab::Backend::native(residual);
// create snippet matrices
std::cout << "create sparse matrices for snippets and subdomains" << std::endl;
auto pvolume_snippet_matrices = std::shared_ptr<std::vector<std::shared_ptr<ISTLM>>>( new std::vector<std::shared_ptr<ISTLM>>(set_up_local_matrices(istlA,pdddm->volumesnippetnumber_local_to_global,pdddm->volumesnippetnumber_global_to_local)) );
for (auto matptr : *pvolume_snippet_matrices) *matptr = 0.0; // clear subdomain matrices
std::shared_ptr<std::vector<std::shared_ptr<ISTLM>>> pskeleton_snippet_matrices; // the empty pointer
pskeleton_snippet_matrices = std::shared_ptr<std::vector<std::shared_ptr<ISTLM>>>( new std::vector<std::shared_ptr<ISTLM>>(set_up_local_matrices(istlA,pdddm->skeletonsnippetnumber_local_to_global,pdddm->skeletonsnippetnumber_global_to_local)) );
for (auto matptr : *pskeleton_snippet_matrices) *matptr = 0.0; // clear subdomain matrices
// make wrapped local operator
std::cout << "assembling subdomain matrices" << std::endl;
using WLOP = SnippetAssemblerWrapper<LocalOperator,ISTLM,VectorGridFunctionSpace>;
WLOP wlop(vgfs,lop,pdddm,pvolume_snippet_matrices,pskeleton_snippet_matrices);
using WGOP = Dune::PDELab::GridOperator<VectorGridFunctionSpace,VectorGridFunctionSpace,WLOP,
Dune::PDELab::ISTL::BCRSMatrixBackend<>,double,double,double,
ConstraintsContainer,ConstraintsContainer>;
WGOP wgop(vgfs,cc,vgfs,cc,wlop,mb);
Dune::Timer timer;
timer.reset();
wgop.jacobian(displacement,A); // assembles the global stiffness matrix and the snippets
auto assemble_time = timer.elapsed();
std::cout << "SNIPPET ASSEMBLE TIME=" << assemble_time << std::endl;
// global vector with dirichlet boundary conditions
// contains 0.0 for dof on Dirichlet boundary, otherwise 1.0
Vector dirichlet_mask(vgfs);
dirichlet_mask = 1.0;
set_constrained_dofs(cc,0.0,dirichlet_mask); // dirichlet dofs are zero, others are 1, so we can use it as a mask!
std::shared_ptr<ISTLV> pistldirichlet_mask = stackobject_to_shared_ptr(Dune::PDELab::Backend::native(dirichlet_mask));
// symmetric elimination of Dirichlet boundary conditions; so we can check symmetry
auto blocksize = (*pistldirichlet_mask)[0].size();
for (size_t i=0; i<istlA.N(); i++)
for (size_t compi=0; compi<blocksize; compi++)
if ((*pistldirichlet_mask)[i][compi]==1.0) // this dof is NOT on the Dirichlet boundary
{
// non dirchlet row; eliminate Dirichlet columns
auto cIt = istlA[i].begin();
auto cEndIt = istlA[i].end();
for (; cIt!=cEndIt; ++cIt)
for (size_t compj=0; compj<blocksize; compj++)
if ((*pistldirichlet_mask)[cIt.index()][compj]==0.0)
(*cIt)[compi][compj] = 0.0;
}
bool sym = is_symmetric(istlA,"A",true,1e-12);
std::cout << "matrix is symmetric: " << sym << std::endl;
// set up ISTL solver
int itmax = 100000;
using Operator = Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV>;
Operator op(istlA);
// Dune::SeqILU<ISTLM,ISTLV,ISTLV> prec(istlA,0,1.0,false);
Dune::Richardson<ISTLV,ISTLV> prec(1.0);
Dune::CGSolver<ISTLV> solver(op,prec,1e-8,itmax,2);
//Dune::RestartedGMResSolver<ISTLV> solver(op,prec,1e-8,50,itmax,2);
Dune::InverseOperatorResult stat;
Vector v(vgfs);
ISTLV& istlv = Dune::PDELab::Backend::native(v);
v = 1.0;
set_constrained_dofs(cc,0.0,v); // dirichlet dofs are zero, others are 1, so we can use it as a mask!
solver.apply(istlv,istlr,stat);
istldisplacement -= istlv;
//==================
// End Geneo stuff
//==================
// The old solver
// slp.apply();
// Visualization
Dune::SubsamplingVTKWriter vtkwriter(gv, Dune::RefinementIntervals(2));
// Add the displacement field visualization
Dune::PDELab::addSolutionToVTKWriter(vtkwriter, vgfs, displacement , Dune::PDELab::vtk::DefaultFunctionNameGenerator("Displacement"));
// Add von Mises stress visualization
using P1FEM = Dune::PDELab::PkLocalFiniteElementMap<EntitySet, double, double, 1>;
using P1GFS = Dune::PDELab::GridFunctionSpace<EntitySet, P1FEM, Dune::PDELab::NoConstraints, LeafVectorBackend>;
using StressVector = Dune::PDELab::Backend::Vector<P1GFS, double>;
P1FEM p1fem(es);
P1GFS p1gfs(es, p1fem);
StressVector stress_container(p1gfs);
VonMisesStressGridFunction<Vector, 2> stress(displacement, material);
Dune::PDELab::interpolate(stress, p1gfs, stress_container);
Dune::PDELab::addSolutionToVTKWriter(vtkwriter, p1gfs, stress_container, Dune::PDELab::vtk::DefaultFunctionNameGenerator("von-Mises stress"));
// Write visualization to file
vtkwriter.write("output", Dune::VTK::appendedraw);
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment