Skip to content
Snippets Groups Projects
Commit 7e310e6f authored by Dominic Kempf's avatar Dominic Kempf
Browse files

Finish transition to BlockLab

parent f33ba658
No related branches found
No related tags found
No related merge requests found
Showing
with 97 additions and 795 deletions
......@@ -28,11 +28,10 @@ link_libraries(stdc++fs)
add_subdirectory(doc)
add_subdirectory(operators)
add_subdirectory(apps)
add_subdirectory(app)
add_subdirectory(cmake/modules)
add_subdirectory(dune)
add_subdirectory(python)
add_subdirectory(test)
add_subdirectory(experiments)
# finalize the dune project, e.g. generating config.h etc.
......
add_executable(structures EXCLUDE_FROM_ALL structures.cc)
#include"config.h"
/** This is the main app that allows to access all of BlockLab's core
* capabilities through one executable. Of course, it takes horribly
* long to compile. For downstream projects, it makes much more sense
* to provide their own app, that follows the implementation patterns
* used in this app, but only includes the necessary grid provider,
* vector providers and blocks.
*/
#include<dune/blocklab/construction/combinatorics.hh>
#include<dune/blocklab/grids.hh>
#include<dune/blocklab/utilities/tuple.hh>
#include<dune/blocklab/vectors.hh>
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertree.hh>
#include<dune/common/parametertreeparser.hh>
#include<dune/pdelab/common/partitionviewentityset.hh>
#include<dune/structures/elasticity.hh>
#include<dune/structures/gridprovider.hh>
#include<dune/structures/material.hh>
#include<dune/structures/onetoone.hh>
#include<dune/structures/visualization.hh>
#include<dune/structures/reinforcedoperator.hh>
#include<memory>
#include<tuple>
template<typename... P>
struct ParameterGatherer
{
using Materials = std::tuple<std::shared_ptr<ElasticMaterialBase<Dune::PDELab::OverlappingEntitySet<typename P::Grid::Traits::LeafGridView>, double>>...>;
};
int main(int argc, char** argv)
{
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
// Parse the ini file
Dune::ParameterTree config;
Dune::ParameterTreeParser::readOptions(argc, argv, config);
Dune::ParameterTreeParser::readINITree(argv[1], config, false);
auto grid_providers = std::make_tuple(
// std::make_pair("cell3d", [](const auto& c){ return std::make_shared<StructuresGridProvider<3>>(c); })//,
std::make_pair("structuredug2d", [](const auto& c){ return std::make_shared<Dune::BlockLab::StructuredSimplexGridProvider<Dune::UGGrid<2>>>(c); })//,
// std::make_pair("structuredug3d", [](const auto& c){ return std::make_shared<Dune::BlockLab::StructuredSimplexGridProvider<Dune::UGGrid<3>>>(c); })
);
auto vector_providers = std::make_tuple(
std::make_pair("p1fem",
[](const auto& c, auto gp){
using GridProvider = typename decltype(gp)::element_type;
auto leaf = std::make_shared<Dune::BlockLab::PkFemVectorProvider<GridProvider, 1>>(gp);
return Dune::BlockLab::fieldProvider<GridProvider::Grid::dimension>(leaf);
})//,
// std::make_pair("p2fem",
// [](const auto& c, auto gp){
// using GridProvider = typename decltype(gp)::element_type;
// auto leaf = std::make_shared<Dune::BlockLab::PkFemVectorProvider<GridProvider, 2>>(gp);
// return Dune::BlockLab::powerProvider<GridProvider::Grid::dimension>(leaf);
// })
);
// The registration function that we are using
auto reg = [](auto& ctx){
Dune::BlockLab::registerBuiltinBlocks(ctx);
// Add the structures-specific blocks to the solver
ctx.template registerBlock<ElasticityOperatorBlock>("elasticity_operator");
ctx.template registerBlock<MaterialInitializationBlock>("material");
ctx.template registerBlock<FibreReinforcedElasticityOperatorBlock>("reinforced_operator");
ctx.template registerBlock<ElasticityMassOperatorBlock>("mass_operator");
ctx.template registerBlock<OneToOneMappingCheckerBlock>("onetoone");
ctx.template registerBlock<FibreDistanceVisualizationBlock>("vis_fibredistance");
ctx.template registerBlock<PhysicalEntityVisualizationBlock>("vis_physical");
ctx.template registerBlock<VonMisesStressVisualizationBlock>("vis_vonmises");
};
// Construct additional types to inject into the parameter system
using MaterialTuple = ParameterGatherer<
StructuresGridProvider<3>,
Dune::BlockLab::StructuredSimplexGridProvider<Dune::UGGrid<2>>,
Dune::BlockLab::StructuredSimplexGridProvider<Dune::UGGrid<3>>
>::Materials;
using ParameterTuple = Dune::BlockLab::tuple_cat_t<std::tuple<std::shared_ptr<std::vector<int>>>, MaterialTuple>;
Dune::BlockLab::instantiate_combinatorics<3, ParameterTuple>(helper, config, grid_providers, vector_providers, reg);
return 0;
}
#add_subdirectory(dynamics)
add_subdirectory(universal)
add_executable(dynamicsapp dynamicsapp.cc)
add_dependencies(dynamicsapp operators)
dune_symlink_to_source_files(FILES linear_dyncube.ini
neohookean_dyncube.ini
stvenant_dyncube.ini
)
#include"config.h"
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertree.hh>
#include<dune/structures/elasticity.hh>
#include<dune/structures/gridconstruction.hh>
#include<dune/structures/solverconstruction.hh>
// This is necessary to enable C++ string literals like "foo"s
// We need these in some places.
using namespace std::literals;
int main(int argc, char** argv)
{
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
// Parse the ini file
Dune::ParameterTree params;
Dune::ParameterTreeParser::readINITree(argv[1], params);
auto [grid, es, physical] = construct_grid(helper, params.sub("grid"), argv);
auto [x, cc] = elastodynamics_setup(es);
using V = std::remove_reference<decltype(*x)>::type;
ElastodynamicsConstructionContext<V> ctx(helper, params, es, physical);
auto solver = ctx.construct(params.sub("solver"));
solver->apply(x, cc);
return 0;
}
[solver]
steps = parameter, material, interpolation, constraints1, timeloop1, constraints2, timeloop2
[grid]
type = structured
N = 10 10 10
[newton]
VerbosityLevel = 3
Reduction = 1e-12
MaxIterations = 100
[parameter]
value = 0.2
datatype = double
name = shear
[material]
materials = cube
[material.cube]
group = 0
youngs_modulus = 10000
poisson_ratio = 0.4
model = linear
[interpolation]
functions = 0.0
[constraints1]
type = constraints
functions = (z < 1e-08) || (z > 1.0 - 1e-8)
[timeloop1]
type = timeloop
Tstart = 0.0
Tend = 1.0
dt = 0.2
steps = elastodynamics_bc
[elastodynamics_bc]
type = elastodynamics
functions = shear * z * time, 0.0, 0.0, shear * z, 0.0, 0.0
[constraints2]
type = constraints
functions = z < 1e-08
[timeloop2]
type = timeloop
Tstart = 1.0
Tend = 1.3
dt = 0.001
steps = elastodynamics, visualization
[visualization]
name = shear_linear
path = linvtu
steps = vis_solution
\ No newline at end of file
[solver]
steps = parameter, material, interpolation, constraints1, timeloop1, constraints2, timeloop2
[grid]
type = structured
N = 10 10 10
[newton]
VerbosityLevel = 3
Reduction = 1e-12
MaxIterations = 100
[parameter]
value = 0.2
datatype = double
name = shear
[material]
materials = cube
[material.cube]
group = 0
youngs_modulus = 10000
poisson_ratio = 0.4
model = neohookean
[interpolation]
functions = 0.0
[constraints1]
type = constraints
functions = (z < 1e-08) || (z > 1.0 - 1e-8)
[timeloop1]
type = timeloop
Tstart = 0.0
Tend = 1.0
dt = 0.2
steps = elastodynamics_bc
[elastodynamics_bc]
type = elastodynamics
functions = shear * z * time, 0.0, 0.0, shear * z, 0.0, 0.0
[constraints2]
type = constraints
functions = z < 1e-08
[timeloop2]
type = timeloop
Tstart = 1.0
Tend = 1.3
dt = 0.001
steps = elastodynamics, visualization
[visualization]
name = shear_neohookean
path = neovtu
steps = vis_solution
[solver]
steps = parameter, material, interpolation, constraints1, timeloop1, constraints2, timeloop2
[grid]
type = structured
N = 10 10 10
[newton]
VerbosityLevel = 3
Reduction = 1e-12
MaxIterations = 100
[parameter]
value = 0.2
datatype = double
name = shear
[material]
materials = cube
[material.cube]
group = 0
youngs_modulus = 10000
poisson_ratio = 0.4
model = stvenantkirchhoff
[interpolation]
functions = 0.0
[constraints1]
type = constraints
functions = (z < 1e-08) || (z > 1.0 - 1e-8)
[timeloop1]
type = timeloop
Tstart = 0.0
Tend = 1.0
dt = 0.2
steps = elastodynamics_bc
[elastodynamics_bc]
type = elastodynamics
functions = shear * z * time, 0.0, 0.0, shear * z, 0.0, 0.0
[constraints2]
type = constraints
functions = z < 1e-08
[timeloop2]
type = timeloop
Tstart = 1.0
Tend = 1.3
dt = 0.001
steps = elastodynamics, visualization
[visualization]
name = stvenantkirchhoff
path = svkvtu
steps = vis_solution
add_executable(universalapp universalapp.cc)
add_dependencies(universalapp operators)
dune_symlink_to_source_files(FILES fibrecube.ini
iterative_cube_compress.ini
linearized_cube_compress.ini
quasistatic.ini
robert.ini
robert2.ini
simplecell.ini
twod_beam.ini
)
[solver]
steps = material, interpolation, constraints, elasticity_operator, newton, visualization, onetoone
[grid]
type = structures
filename = fibrecube.msh
[grid.cytoplasm]
shape = box
physical = 0
[grid.fibres]
fibres = one
[grid.fibres.one]
shape = overnucleus
radius = 0.02
start = -0.5 0.0 -0.5
middle = 0.0 0.0 0.1
end = 0.5 0.0 -0.5
slope = 0.5
meshwidth = 0.01
physical = 2
[grid.export.geo]
enabled = 1
[grid.export.vtk]
enabled = 1
[material]
materials = cyto, fibre
[material.cyto]
group = 0
youngs_modulus = 1000
poisson_ratio = 0.4
model = linear
[material.cyto.prestress]
type = none
#isotropic
scale = 1000
[material.fibre]
group = 2
youngs_modulus = 5000
poisson_ratio = 0.4
model = linear
[material.fibre.prestress]
type = curved_fibre
fibre = grid.fibres.one
scale = 5000
[interpolation]
functions = 0.0
[constraints]
functions = z < -0.5 + 1e-08
[newton]
operator = elasticity_operator
[visualization]
name = fibrecube
steps = vis_solution, vis_vonmises, vis_fibredistance
[vis_fibredistance]
key = material.fibre.prestress
\ No newline at end of file
[solver]
steps = parameter, material, elasticity_operator, interpolation, constraints, compression, visualization
[grid]
type = structured
N = 8 8 8
[newton]
VerbosityLevel = 3
Reduction = 1e-12
MaxIterations = 100
[parameter]
value = 0.5
datatype = double
name = maxcompression
[material]
materials = cube
[material.cube]
group = 0
youngs_modulus = 10000
poisson_ratio = 0.4
model = neohookean
[constraints]
functions = (z < 1e-08) || (z > 1.0 - 1e-8)
[interpolation]
functions = 0.0
[compression]
type = continuousvariation
name = compression
iterations = 10
start = 0.0
end = 1.0
steps = transformation, newton, onetoone
[newton]
operator = elasticity_operator
[transformation]
functions = ux, uy, - compression * (1.0 - maxcompression) * z
[visualization]
name = iterative_compress
steps = vis_solution, vis_vonmises
[solver]
steps = parameter, material, elasticity_operator, interpolation, constraints, linearization, visualization
[grid]
type = structured
N = 8 8 8
[newton]
VerbosityLevel = 3
Reduction = 1e-12
MaxIterations = 100
[parameter]
value = 0.7
datatype = double
name = maxcompression
[material]
materials = cube
[material.cube]
group = 0
youngs_modulus = 10000
poisson_ratio = 0.4
model = stvenantkirchhoff
[constraints]
functions = (z < 1e-08) || (z > 1.0 - 1e-8)
[interpolation]
functions = 0.0, 0.0, -(1.0 - maxcompression) * z
[linearization]
type = discretematerialvariation
name = model
key = cube.model
values = linear, stvenantkirchhoff
steps = newton, onetoone
[newton]
operator = elasticity_operator
[visualization]
name = linearized_compress
steps = vis_solution, vis_vonmises
[solver]
steps = material, elasticity_operator, elasticity_mass_operator, interpolation, constraints, timeloop
[grid]
type = structured
N = 8 8 8
[material]
materials = cube
[material.cube]
group = 0
youngs_modulus = 10000
poisson_ratio = 0.4
model = linear
[constraints]
functions = z < 1e-08
[interpolation]
functions = 0.2 * z, 0.0, 0.0
[timeloop]
Tstart = 0.0
Tend = 0.1
dt = 0.001
steps = onestep, visualization
[onestep]
spatial_operator = elasticity_operator
temporal_operator = elasticity_mass_operator
[visualization]
instationary = 1
name = quasistatic
path = quasistatic
steps = vis_solution
[solver]
steps = material, interpolation, constraints, elasticity_operator, newton, visualization
[grid]
type = cell
filename = robert
scaling = 1e-6
[grid.cytoplasm]
meshwidth = 0.5
physical = 0
shape = box
lowerleft = 0.0 0.0 0.0
size = 10.0 10.0 5.0
[grid.fibres]
fibres = first
[grid.fibres.first]
shape = cylinder
meshwidth = 0.05
start = 0.0 5.0 2.5
end = 10.0 5.0 2.5
radius = 0.1
physical = 1
[grid.export.vtk]
enabled = 1
[material]
materials = cyto, fibre
[material.cyto]
group = 0
youngs_modulus = 1000
poisson_ratio = 0.4
model = linear
#[material.cyto.prestress]
#type = isotropic
#scale = 500
[material.fibre]
group = 1
youngs_modulus = 100000
poisson_ratio = 0.4
model = linear
#[material.fibre.prestress]
#type = straight_fibre
#direction = 1.0 0.0 0.0
#scale = 1000
[interpolation]
functions = 0.0
[constraints]
functions = z < 1e-10
[elasticity_operator]
traction = (group == 1 ? 1.0 : 0.0) * ((x < 1e-8) ? -1.0 : 1.0) * 100000, 0.0, 0.0
[newton]
operator = elasticity_operator
[visualization]
name = robert
steps = vis_solution, vis_physicalentity, vis_vonmises
[solver]
steps = material, elasticity_operator, interpolation, constraints, newton, visualization, onetoone
[grid]
type = cell
filename = robert2
scaling = 1e-6
[grid.cytoplasm]
meshwidth = 0.5
physical = 0
shape = box
lowerleft = 0.0 0.0 0.0
size = 10.0 10.0 2.0
[grid.fibres]
fibres = first, second
[grid.fibres.first]
shape = overnucleus
meshwidth = 0.1
start = 0.07 0.07 0.0
middle = 5.0 5.0 1.5
end = 9.93 9.93 0.0
slope = 0.5
radius = 0.1
physical = 1
[grid.fibres.second]
shape = overnucleus
meshwidth = 0.1
start = 0.07 9.93 0.0
middle = 5.0 5.0 1.25
end = 9.93 0.07 0.0
slope = 0.5
radius = 0.1
physical = 2
[grid.export.vtk]
enabled = 1
[material]
materials = cyto, fibre1, fibre2
[material.cyto]
group = 0
youngs_modulus = 1000
poisson_ratio = 0.45
model = linear
[material.cyto.prestress]
type = isotropic
scale = 500
[material.fibre1]
group = 1
youngs_modulus = 300000
poisson_ratio = 0.45
model = linear
[material.fibre1.prestress]
type = curved_fibre
scale = 300000
fibre = grid.fibres.first
[material.fibre2]
group = 2
youngs_modulus = 300000
poisson_ratio = 0.45
model = linear
[material.fibre2.prestress]
type = curved_fibre
scale = 300000
fibre = grid.fibres.second
[interpolation]
functions = 0.0
[constraints]
functions = group > 0, group > 0, z<1e-12
[newton]
operator = elasticity_operator
[visualization]
steps = vis_solution, vis_vonmises, vis_physicalentity
name = robert2
[solver]
steps = material, interpolation, constraints, elasticity_operator, newton, visualization, onetoone
[grid]
type = cell
filename = spread.msh
scaling = 1e-6
[grid.cytoplasm]
shape = spread
radius = 12.5
height = 7.5
slope = 0.5
meshwidth = 0.5
physical = 0
[grid.nucleus]
shape = ellipsoid
center = 0.0 0.0 3.75
radii = 4.5 4.5 2.25
physical = 1
meshwidth = 0.5
[grid.fibres]
fibres = one
[grid.fibres.one]
shape = overnucleus
radius = 0.1
start = -12.5 0.0 -0.25
middle = 0.0 0.0 6.15
end = 12.5 0.0 -0.25
slope = 0.5
meshwidth = 0.05
physical = 2
[grid.export.geo]
enabled = 1
[grid.export.vtk]
enabled = 1
[material]
materials = cyto, nucleus, fibre
[material.cyto]
group = 0
youngs_modulus = 1000
poisson_ratio = 0.4
model = linear
[material.cyto.prestress]
type = isotropic
scale = 500
[material.nucleus]
group = 1
youngs_modulus = 10000
poisson_ratio = 0.4
model = linear
[material.fibre]
group = 2
youngs_modulus = 5000
poisson_ratio = 0.4
model = linear
[material.fibre.prestress]
type = curved_fibre
fibre = grid.fibres.one
scale = 5000
[interpolation]
functions = 0.0
[constraints]
functions = z < 1e-08
[elasticity_operator]
[newton]
operator = elasticity_operator
[probe]
position = 0.0 0.0 1.5e-6
name = nucleus_bottom
displace = 1
[visualization]
name = simplecell
steps = vis_solution, vis_vonmises, vis_physicalentity
[solver]
steps = material, interpolation, constraints, elasticity_operator, fibre_operator, linearsolver, visualization
degree = 2
[grid]
dimension = 2
type = structured
lowerleft = 0.0 0.0
upperright = 4.0 1.0
N = 40 10
[grid.fibres]
fibres = first
[grid.fibres.first]
shape = cylinder
start = -0.5 0.501
end = 4.5 0.501
youngs_modulus = 1e6
radius = 0.1
[material]
materials = bulk
[material.bulk]
model = linear
youngs_modulus = 300
poisson_ratio = 0.3333
[interpolation]
functions = 0.0
[constraints]
functions = x < 1e-8
[fibre_operator]
force = 0, -1
#force = 100, 0
[elasticity_operator]
force = 0, -1
[linearsolver]
operator = elasticity_operator
#operator = fibre_operator
[visualization]
steps = vis_solution
name = twodbeam
#include"config.h"
#include<dune/common/parallel/mpihelper.hh>
#include<dune/common/parametertree.hh>
#include<dune/structures/elasticity.hh>
#include<dune/structures/gridconstruction.hh>
#include<dune/structures/gridreorder.hh>
#include<dune/structures/solverconstruction.hh>
#include<memory>
#include<tuple>
template<int dim, int degree>
void apply(Dune::MPIHelper& helper, const Dune::ParameterTree& params, char** argv)
{
auto [grid, es, physical] = construct_grid<dim>(helper, params.sub("grid"), argv);
std::tie(grid, es, physical) = reorder_grid(grid, es, physical, params.sub("grid"));
auto [x, cc, force, force_cc, traction, traction_cc] = elasticity_setup<degree>(es);
using V = typename std::remove_reference<decltype(*x)>::type;
using V_F = typename std::remove_reference<decltype(*force)>::type;
using V_T = typename std::remove_reference<decltype(*traction)>::type;
ConstructionContext<V, V_F, V_T> ctx(helper, params, grid, es, physical);
ctx.setVectors(x, force, traction);
ctx.setConstraintsContainers(cc, force_cc, traction_cc);
auto solver = ctx.construct(params.sub("solver"));
solver->apply();
}
int main(int argc, char** argv)
{
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
// Parse the ini file
Dune::ParameterTree params;
Dune::ParameterTreeParser::readOptions(argc, argv, params);
Dune::ParameterTreeParser::readINITree(argv[1], params, false);
// Dispatch on grid dimension
auto dim = params.get<int>("grid.dimension", 3);
auto degree = params.get<int>("solver.degree", 1);
if (dim == 2)
{
if (degree == 1)
apply<2, 1>(helper, params, argv);
if (degree == 2)
apply<2, 2>(helper, params, argv);
}
if (dim == 3)
{
if (degree == 1)
apply<3, 1>(helper, params, argv);
if (degree == 2)
apply<3, 2>(helper, params, argv);
}
return 0;
}
find_package(muparser REQUIRED)
dune_register_package_flags(LIBRARIES muparser::muparser)
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment