diff --git a/.github/workflows/building.yml b/.github/workflows/building.yml index 166e60d..9af3f96 100644 --- a/.github/workflows/building.yml +++ b/.github/workflows/building.yml @@ -22,10 +22,6 @@ jobs: cd dealii-adapter && \ git fetch origin ${{ github.ref }} && \ git checkout FETCH_HEAD && \ - cd linear_elasticity && \ - cmake . && \ - make && \ - cd ../nonlinear_elasticity && \ cmake . && \ make"; diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..0cbbdbd --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,76 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +FIND_PACKAGE(deal.II 9.1 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +# Enable a switchable dimension choice +IF (NOT DEFINED DIM) + SET(DIM 2) +ENDIF() +ADD_DEFINITIONS(-DDIM=${DIM}) + +# Set the target and the target source +SET( TARGET "elasticity" ) +SET( TARGET_SRC ${TARGET}.cc include/adapter/parameters.cc) + +DEAL_II_INITIALIZE_CACHED_VARIABLES() + +IF(NOT CMAKE_BUILD_TYPE) + SET(CMAKE_BUILD_TYPE "RELEASE") + MESSAGE(STATUS "No build type specified. Building in ${CMAKE_BUILD_TYPE} mode.") +ENDIF() + +# Set the include directory and the name of the project +PROJECT(${TARGET}) +INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/include) + +# Set the executable +ADD_EXECUTABLE(${TARGET} ${TARGET_SRC}) + +# Use deal.II macros for setting up the target +DEAL_II_SETUP_TARGET(${TARGET}) +DEAL_II_INITIALIZE_CACHED_VARIABLES() + +# Query the git information and set it in the source +DEAL_II_QUERY_GIT_INFORMATION() +SET_PROPERTY(TARGET ${TARGET} APPEND PROPERTY COMPILE_DEFINITIONS + GIT_BRANCH="${GIT_BRANCH}" + GIT_REVISION="${GIT_REVISION}" + GIT_SHORTREV="${GIT_SHORTREV}") + + +FIND_PACKAGE(precice REQUIRED + HINTS ${PRECICE_DIR} $ENV{precice_DIR}) +TARGET_LINK_LIBRARIES(${TARGET} precice::precice) + +MESSAGE(STATUS "Using the preCICE version found at ${precice_CONFIG}") + +# add the individual solver +ADD_SUBDIRECTORY(source) + +# ...and link them +TARGET_LINK_LIBRARIES(${TARGET} nonlinear_elasticity linear_elasticity) + + +# +# Custom "debug" and "release" make targets: +# +ADD_CUSTOM_TARGET(debug +COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Debug ${CMAKE_SOURCE_DIR} +COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all +COMMENT "Switch CMAKE_BUILD_TYPE to Debug" +) + +ADD_CUSTOM_TARGET(release +COMMAND ${CMAKE_COMMAND} -DCMAKE_BUILD_TYPE=Release ${CMAKE_SOURCE_DIR} +COMMAND ${CMAKE_COMMAND} --build ${CMAKE_BINARY_DIR} --target all +COMMENT "Switch CMAKE_BUILD_TYPE to Release" +) diff --git a/README.md b/README.md index 3b197e6..f36f035 100644 --- a/README.md +++ b/README.md @@ -6,10 +6,10 @@ Coupled structural solvers written with the C++ finite element library deal.II: -- `linear_elasticity` contains a linear-elastic solver based on the [step-8 tutorial program](https://www.dealii.org/developer/doxygen/deal.II/step_8.html) of deal.II -- `nonlinear_elasticity` contains a nonlinear elastic solver, which builds on previous work of Jean-Paul Pelteret and Andrew McBride in their deal.II code gallery program '[Quasi-Static Finite-Strain Compressible Elasticity](https://www.dealii.org/developer/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html).' This solver supports shared-memory parallelization. +- `source/linear_elasticity` contains a linear-elastic solver based on the [step-8 tutorial program](https://www.dealii.org/developer/doxygen/deal.II/step_8.html) of deal.II +- `source/nonlinear_elasticity` contains a nonlinear elastic solver, which builds on previous work of Jean-Paul Pelteret and Andrew McBride in their deal.II code gallery program '[Quasi-Static Finite-Strain Compressible Elasticity](https://www.dealii.org/developer/doxygen/deal.II/code_gallery_Quasi_static_Finite_strain_Compressible_Elasticity.html).' This solver supports shared-memory parallelization. -Applied coupling functionalities have been separated and can be found in the `adapter` directory. +Applied coupling functionalities have been separated and can be found in the `include/adapter` directory. ## Start here Our [wiki](https://www.precice.org/adapter-dealii-overview.html) will help you start. If you are missing something, [let us know](https://www.precice.org/resources/#contact). diff --git a/elasticity.cc b/elasticity.cc new file mode 100644 index 0000000..006cdef --- /dev/null +++ b/elasticity.cc @@ -0,0 +1,101 @@ +#include + +#include "source/linear_elasticity/include/linear_elasticity.h" +#include "source/nonlinear_elasticity/include/nonlinear_elasticity.h" + +int +main(int argc, char **argv) +{ + using namespace dealii; + +#ifdef DEAL_II_WITH_MPI + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); +#endif + + try + { + deallog.depth_console(0); + static const unsigned int n_threads = MultithreadInfo::n_threads(); + + // Query adapter and deal.II info + const std::string adapter_info = + GIT_SHORTREV == std::string("") ? + "unknown" : + (GIT_SHORTREV + std::string(" on branch ") + GIT_BRANCH); + const std::string dealii_info = + DEAL_II_GIT_SHORTREV == std::string("") ? + "unknown" : + (DEAL_II_GIT_SHORTREV + std::string(" on branch ") + + DEAL_II_GIT_BRANCH); + + std::cout + << "-----------------------------------------------------------------------------" + << std::endl + << "-- . running with " << n_threads << " thread" + << (n_threads == 1 ? "" : "s") << std::endl; + + std::cout << "-- . adapter revision " << adapter_info << std::endl; + std::cout << "-- . deal.II " << DEAL_II_PACKAGE_VERSION + << " (revision " << dealii_info << ")" << std::endl; + std::cout + << "-----------------------------------------------------------------------------" + << std::endl + << std::endl; + + // Store the name of the parameter file + const std::string parameter_file = argc > 1 ? argv[1] : "parameters.prm"; + + // Extract case path for the output directory + size_t pos = parameter_file.find_last_of("/"); + std::string case_path = + std::string::npos == pos ? "" : parameter_file.substr(0, pos + 1); + + // Query solver type from the parameter file + ParameterHandler prm; + Parameters::Solver solver; + solver.add_output_parameters(prm); + prm.parse_input(parameter_file, "", true); + + if (solver.model == "neo-Hookean") // nonlinear + { + Nonlinear_Elasticity::Solid solid(case_path, parameter_file); + solid.run(); + } + else if (solver.model == "linear") // linear + { + Linear_Elasticity::ElastoDynamics elastic_solver(case_path, + parameter_file); + elastic_solver.run(); + } + else + AssertThrow(false, ExcNotImplemented()) + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/adapter/adapter.h b/include/adapter/adapter.h similarity index 99% rename from adapter/adapter.h rename to include/adapter/adapter.h index ace85cc..508cc9a 100644 --- a/adapter/adapter.h +++ b/include/adapter/adapter.h @@ -9,10 +9,10 @@ #include #include +#include +#include #include -#include "dof_tools_extension.h" -#include "time.h" namespace Adapter { diff --git a/adapter/dof_tools_extension.h b/include/adapter/dof_tools_extension.h similarity index 100% rename from adapter/dof_tools_extension.h rename to include/adapter/dof_tools_extension.h diff --git a/include/adapter/parameters.cc b/include/adapter/parameters.cc new file mode 100644 index 0000000..dc73797 --- /dev/null +++ b/include/adapter/parameters.cc @@ -0,0 +1,202 @@ +#include + +namespace Parameters +{ + void + Time::add_output_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + prm.add_parameter("End time", end_time, "End time", Patterns::Double()); + + prm.add_parameter("Time step size", + delta_t, + "Time step size", + Patterns::Double()); + + prm.add_parameter("Output interval", + output_interval, + "Write results every x timesteps", + Patterns::Integer(0)); + } + prm.leave_subsection(); + } + + + void + System::add_output_parameters(ParameterHandler &prm) + { + prm.enter_subsection("System properties"); + { + prm.add_parameter("Shear modulus", + mu, + "Shear modulus", + Patterns::Double()); + + prm.add_parameter("Poisson's ratio", + nu, + "Poisson's ratio", + Patterns::Double(-1.0, 0.5)); + + prm.add_parameter("rho", rho, "density", Patterns::Double(0.0)); + + prm.add_parameter("body forces", + body_force, + "body forces x,y,z", + Patterns::List(Patterns::Double())); + } + prm.leave_subsection(); + + lambda = 2 * mu * nu / (1 - 2 * nu); + } + + + void + Solver::add_output_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Solver"); + { + prm.add_parameter("Model", + model, + "Structural model to be used: linear or neo-Hookean", + Patterns::Selection("linear|neo-Hookean")); + + prm.add_parameter("Solver type", + type_lin, + "Linear solver: CG or Direct", + Patterns::Selection("CG|Direct")); + + prm.add_parameter( + "Residual", + tol_lin, + "CG solver residual (multiplied by residual norm, ignored if Model == linear)", + Patterns::Double(0.0)); + + prm.add_parameter( + "Max iteration multiplier", + max_iterations_lin, + "Max CG solver iterations (multiples of the system matrix size)", + Patterns::Double(0.0)); + + prm.add_parameter( + "Max iterations Newton-Raphson", + max_iterations_NR, + "Number of Newton-Raphson iterations allowed (ignored if Model == linear)", + Patterns::Integer(0)); + + prm.add_parameter( + "Tolerance force", + tol_f, + "Force residual tolerance for non-linear iteration (ignored if Model == linear)", + Patterns::Double(0.0)); + + prm.add_parameter( + "Tolerance displacement", + tol_u, + "Displacement error tolerance for non-linear iteration (ignored if Model == linear)", + Patterns::Double(0.0)); + } + prm.leave_subsection(); + } + + + void + Discretization::add_output_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Discretization"); + { + prm.add_parameter("Polynomial degree", + poly_degree, + "Polynomial degree of the FE system", + Patterns::Integer(0)); + + prm.add_parameter("theta", + theta, + "Time integration scheme", + Patterns::Double(0, 1)); + + prm.add_parameter("beta", beta, "Newmark beta", Patterns::Double(0, 0.5)); + + prm.add_parameter("gamma", + gamma, + "Newmark gamma", + Patterns::Double(0, 1)); + } + prm.leave_subsection(); + } + + + void + PreciceAdapterConfiguration::add_output_parameters(ParameterHandler &prm) + { + prm.enter_subsection("precice configuration"); + { + prm.add_parameter("Scenario", + scenario, + "Cases: FSI3 or PF for perpendicular flap", + Patterns::Selection("FSI3|PF")); + + prm.add_parameter("precice config-file", + config_file, + "Name of the precice configuration file", + Patterns::Anything()); + + prm.add_parameter( + "Participant name", + participant_name, + "Name of the participant in the precice-config.xml file", + Patterns::Anything()); + + prm.add_parameter( + "Mesh name", + mesh_name, + "Name of the coupling mesh in the precice-config.xml file", + Patterns::Anything()); + + prm.add_parameter("Read data name", + read_data_name, + "Name of the read data in the precice-config.xml file", + Patterns::Anything()); + + prm.add_parameter("Write data name", + write_data_name, + "Name of the write data in the precice-config.xml file", + Patterns::Anything()); + + prm.add_parameter("Flap location", + flap_location, + "PF x-location", + Patterns::Double(-3, 3)); + } + prm.leave_subsection(); + } + + + AllParameters::AllParameters(const std::string &input_file) + { + ParameterHandler prm; + + Solver::add_output_parameters(prm); + Discretization::add_output_parameters(prm); + System::add_output_parameters(prm); + Time::add_output_parameters(prm); + PreciceAdapterConfiguration::add_output_parameters(prm); + + prm.parse_input(input_file); + + // Look at the specific type of read data + if ((read_data_name.find("Stress") == 0)) + data_consistent = true; + else if ((read_data_name.find("Force") == 0)) + data_consistent = false; + else + AssertThrow( + false, + ExcMessage( + "Unknown read data type. Please use 'Force' or 'Stress' in the read data naming.")); + + // Optional, if we want to print all parameters in the beginning of the + // simulation + // prm.print_parameters(std::cout,ParameterHandler::Text); + } +} // namespace Parameters diff --git a/include/adapter/parameters.h b/include/adapter/parameters.h new file mode 100644 index 0000000..fc6c6f8 --- /dev/null +++ b/include/adapter/parameters.h @@ -0,0 +1,112 @@ +#ifndef PARAMETERS_H +#define PARAMETERS_H + +#include + +/** + * This class declares all parameters, which can be specified in the parameter + * file. The subsection abut preCICE configurations is directly interlinked + * to the Adapter class. + */ +namespace Parameters +{ + using namespace dealii; + /** + * @brief Time: Specifies simulation time properties + */ + struct Time + { + double end_time = 1; + double delta_t = 0.1; + int output_interval = 1; + + void + add_output_parameters(ParameterHandler &prm); + }; + + /** + * @brief The System struct keeps material properties and body force contributions + */ + struct System + { + double nu = 0.3; + double mu = 1538462; + double lambda = 2307692; + double rho = 1000; + Tensor<1, 3, double> body_force; + + void + add_output_parameters(ParameterHandler &prm); + }; + + + /** + * @brief LinearSolver: Specifies linear solver properties + */ + struct Solver + { + std::string model = "linear"; + std::string type_lin = "Direct"; + double tol_lin = 1e-6; + double max_iterations_lin = 1; + unsigned int max_iterations_NR = 10; + double tol_f = 1e-9; + double tol_u = 1e-6; + + void + add_output_parameters(ParameterHandler &prm); + }; + + + + /** + * @brief Discretization: Specifies parameters for time integration by a + * theta scheme and polynomial degree of the FE system + */ + struct Discretization + { + unsigned int poly_degree = 3; + // For the linear elastic model (theta-scheme) + double theta = 0.5; + // For the nonlinear elastic model (Newmark) + double beta = 0.25; + double gamma = 0.5; + + void + add_output_parameters(ParameterHandler &prm); + }; + + + /** + * @brief PreciceAdapterConfiguration: Specifies preCICE related information. + * A lot of these information need to be consistent with the + * precice-config.xml file. + */ + struct PreciceAdapterConfiguration + { + std::string scenario = "FSI3"; + std::string config_file = "precice-config.xml"; + std::string participant_name = "dealiisolver"; + std::string mesh_name = "dealii-mesh"; + std::string read_data_name = "Stress"; + std::string write_data_name = "Displacement"; + double flap_location = 0.0; + bool data_consistent = true; + + void + add_output_parameters(ParameterHandler &prm); + }; + + + struct AllParameters : public Solver, + public Discretization, + public System, + public Time, + public PreciceAdapterConfiguration + + { + AllParameters(const std::string &input_file); + }; +} // namespace Parameters + +#endif // PARAMETERS_H diff --git a/adapter/time.h b/include/adapter/time_handler.h similarity index 95% rename from adapter/time.h rename to include/adapter/time_handler.h index 80306ba..3e7e486 100644 --- a/adapter/time.h +++ b/include/adapter/time_handler.h @@ -1,5 +1,5 @@ -#ifndef TIME_H -#define TIME_H +#ifndef TIME_HANDLER_H +#define TIME_HANDLER_H #include @@ -81,4 +81,4 @@ namespace Adapter const double delta_t; }; } // namespace Adapter -#endif // TIME_H +#endif // TIME_HANDLER_H diff --git a/linear_elasticity/CMakeLists.txt b/linear_elasticity/CMakeLists.txt deleted file mode 100644 index c7337e6..0000000 --- a/linear_elasticity/CMakeLists.txt +++ /dev/null @@ -1,68 +0,0 @@ -## -# CMake script for the dealii-adapter: -## - -# Set the name of the project and target: -SET(TARGET "linear_elasticity") - -SET(TARGET_SRC - ${TARGET}.cc - ) - -# To make the dimension switch available via cmake -IF (NOT DEFINED DIM) - SET(DIM 2) -ENDIF() -ADD_DEFINITIONS(-DDIM=${DIM}) - -# Change default build type for tutorial cases -IF(NOT CMAKE_BUILD_TYPE) - SET(CMAKE_BUILD_TYPE "Release" CACHE STRING - "Choose the type of build, options are: Debug Release" - FORCE) -ENDIF(NOT CMAKE_BUILD_TYPE) -MESSAGE(STATUS "Build type: ${CMAKE_BUILD_TYPE}") - -# To cleanup result files -SET(CLEAN_UP_FILES - # a custom list of globs, e.g. *.log *.vtk - *.vtk -) -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) - -FIND_PACKAGE(deal.II 9.2.0 QUIET - HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} - ) -IF(NOT ${deal.II_FOUND}) - MESSAGE(FATAL_ERROR "\n" - "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" - "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" - "or set an environment variable \"DEAL_II_DIR\" that contains this path." - ) -ENDIF() - -DEAL_II_INITIALIZE_CACHED_VARIABLES() - -# Store current source directory -SET(TMP ${CMAKE_SOURCE_DIR}) -SET(CMAKE_SOURCE_DIR "..") - -# Query git information -DEAL_II_QUERY_GIT_INFORMATION() - -# Restore source dir -SET(CMAKE_SOURCE_DIR ${TMP}) - -PROJECT(${TARGET} LANGUAGES CXX) - -DEAL_II_INVOKE_AUTOPILOT() - -FIND_PACKAGE(precice REQUIRED) -TARGET_LINK_LIBRARIES(${TARGET} precice::precice) - -SET_PROPERTY(TARGET ${TARGET} APPEND PROPERTY COMPILE_DEFINITIONS - GIT_BRANCH="${GIT_BRANCH}" - GIT_REVISION="${GIT_REVISION}" - GIT_SHORTREV="${GIT_SHORTREV}") diff --git a/linear_elasticity/include/parameter_handling.h b/linear_elasticity/include/parameter_handling.h deleted file mode 100644 index 49e4b8b..0000000 --- a/linear_elasticity/include/parameter_handling.h +++ /dev/null @@ -1,356 +0,0 @@ -#ifndef parameter_handling_h -#define parameter_handling_h - -#include - -namespace Linear_Elasticity -{ - using namespace dealii; - - /** - * This class declares all parameters, which can be specified in the parameter - * file. The subsection abut preCICE configurations is directly interlinked - * to the Adapter class. - */ - namespace Parameters - { - // TODO: Add more parameters - struct Time - { - double delta_t; - double end_time; - int output_interval; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - Time::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - prm.declare_entry("End time", "1", Patterns::Double(), "End time"); - - prm.declare_entry("Time step size", - "0.1", - Patterns::Double(), - "Time step size"); - - prm.declare_entry("Output interval", - "1", - Patterns::Integer(0), - "Write results every x timesteps"); - } - prm.leave_subsection(); - } - - void - Time::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - end_time = prm.get_double("End time"); - delta_t = prm.get_double("Time step size"); - output_interval = prm.get_integer("Output interval"); - } - prm.leave_subsection(); - } - - /** - * @brief Discretization: Specifies parameters for time integration by a - * theta scheme and polynomial degree of the FE system - */ - struct Discretization - { - double theta; - unsigned int poly_degree; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - Discretization::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - prm.declare_entry("theta", - "0.5", - Patterns::Double(0, 1), - "Time integration scheme"); - - prm.declare_entry("Polynomial degree", - "3", - Patterns::Integer(0), - "Polynomial degree of the FE system"); - } - prm.leave_subsection(); - } - - void - Discretization::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - theta = prm.get_double("theta"); - poly_degree = prm.get_integer("Polynomial degree"); - } - prm.leave_subsection(); - } - - /** - * @brief The System struct keeps material properties and body force contributions - */ - struct System - { - double mu; - double lambda; - double rho; - Tensor<1, 3, double> body_force; - - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - System::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - prm.declare_entry("mu", "0.5e6", Patterns::Double(), "mu"); - - prm.declare_entry("lambda", "2e6", Patterns::Double(), "lambda"); - - prm.declare_entry("rho", "1000", Patterns::Double(0.0), "density"); - - prm.declare_entry("body forces", - "0,0,0", - Patterns::List(Patterns::Double()), - "body forces x,y,z"); - } - prm.leave_subsection(); - } - - void - System::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - mu = prm.get_double("mu"); - lambda = prm.get_double("lambda"); - rho = prm.get_double("rho"); - const std::vector body_forces_input = - Utilities::split_string_list(prm.get("body forces")); - for (uint d = 0; d < 3; ++d) - body_force[d] = Utilities::string_to_double(body_forces_input[d]); - } - prm.leave_subsection(); - } - - /** - * @brief LinearSolver: Specifies linear solver properties - */ - struct LinearSolver - { - std::string type_lin; - double tol_lin; - double max_iterations_lin; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - LinearSolver::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - prm.declare_entry("Solver type", - "Direct", - Patterns::Selection("CG|Direct"), - "Linear solver: CG or Direct"); - - prm.declare_entry("Residual", - "1e-6", - Patterns::Double(0.0), - "Linear solver residual (scaled by residual norm)"); - - prm.declare_entry( - "Max iteration multiplier", - "1", - Patterns::Double(0.0), - "Linear solver iterations (multiples of the system matrix size)"); - } - prm.leave_subsection(); - } - - void - LinearSolver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - type_lin = prm.get("Solver type"); - tol_lin = prm.get_double("Residual"); - max_iterations_lin = prm.get_double("Max iteration multiplier"); - } - prm.leave_subsection(); - } - - /** - * @brief PreciceAdapterConfiguration: Specifies preCICE related information. - * A lot of these information need to be consistent with the - * precice-config.xml file. - */ - struct PreciceAdapterConfiguration - { - std::string scenario; - std::string config_file; - std::string participant_name; - std::string mesh_name; - std::string read_data_name; - std::string write_data_name; - double flap_location; - bool data_consistent; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - void - PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - prm.declare_entry("Scenario", - "FSI3", - Patterns::Selection("FSI3|PF"), - "Cases: FSI3 or PF for perpendicular flap"); - prm.declare_entry("precice config-file", - "precice-config.xml", - Patterns::Anything(), - "Name of the precice configuration file"); - prm.declare_entry( - "Participant name", - "dealiisolver", - Patterns::Anything(), - "Name of the participant in the precice-config.xml file"); - prm.declare_entry( - "Mesh name", - "dealii-mesh", - Patterns::Anything(), - "Name of the coupling mesh in the precice-config.xml file"); - prm.declare_entry( - "Read data name", - "received-data", - Patterns::Anything(), - "Name of the read data in the precice-config.xml file"); - prm.declare_entry( - "Write data name", - "calculated-data", - Patterns::Anything(), - "Name of the write data in the precice-config.xml file"); - prm.declare_entry("Flap location", - "0.0", - Patterns::Double(-3, 3), - "PF x-location"); - } - prm.leave_subsection(); - } - - void - PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - scenario = prm.get("Scenario"); - config_file = prm.get("precice config-file"); - participant_name = prm.get("Participant name"); - mesh_name = prm.get("Mesh name"); - read_data_name = prm.get("Read data name"); - write_data_name = prm.get("Write data name"); - flap_location = prm.get_double("Flap location"); - } - prm.leave_subsection(); - // Look at the specific type of read data - if ((read_data_name.find("Stress") == 0)) - data_consistent = true; - else if ((read_data_name.find("Force") == 0)) - data_consistent = false; - else - AssertThrow( - false, - ExcMessage( - "Unknown read data type. Please use 'Force' or 'Stress' in the read data naming.")); - } - - - - struct AllParameters : public LinearSolver, - public Discretization, - public System, - public Time, - public PreciceAdapterConfiguration - - { - AllParameters(const std::string &input_file); - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - AllParameters::AllParameters(const std::string &input_file) - { - ParameterHandler prm; - declare_parameters(prm); - prm.parse_input(input_file); - parse_parameters(prm); - - // Optional, if we want to print all parameters in the beginning of the - // simulation - // prm.print_parameters(std::cout,ParameterHandler::Text); - } - - void - AllParameters::declare_parameters(ParameterHandler &prm) - { - LinearSolver::declare_parameters(prm); - Discretization::declare_parameters(prm); - System::declare_parameters(prm); - Time::declare_parameters(prm); - PreciceAdapterConfiguration::declare_parameters(prm); - } - - void - AllParameters::parse_parameters(ParameterHandler &prm) - { - LinearSolver::parse_parameters(prm); - Discretization::parse_parameters(prm); - System::parse_parameters(prm); - Time::parse_parameters(prm); - PreciceAdapterConfiguration::parse_parameters(prm); - } - } // namespace Parameters -} // namespace Linear_Elasticity - -#endif diff --git a/nonlinear_elasticity/CMakeLists.txt b/nonlinear_elasticity/CMakeLists.txt deleted file mode 100644 index b72a6fc..0000000 --- a/nonlinear_elasticity/CMakeLists.txt +++ /dev/null @@ -1,68 +0,0 @@ -## -# CMake script for the dealii-adapter: -## - -# Set the name of the project and target: -SET(TARGET "nonlinear_elasticity") - -SET(TARGET_SRC - ${TARGET}.cc - ) - -# To make the dimension switch available via cmake -IF (NOT DEFINED DIM) - SET(DIM 2) -ENDIF() -ADD_DEFINITIONS(-DDIM=${DIM}) - -# Change default build type for tutorial cases -IF(NOT CMAKE_BUILD_TYPE) - SET(CMAKE_BUILD_TYPE "Release" CACHE STRING - "Choose the type of build, options are: Debug Release" - FORCE) -ENDIF(NOT CMAKE_BUILD_TYPE) -MESSAGE(STATUS "Build type: ${CMAKE_BUILD_TYPE}") - -# To cleanup result files -SET(CLEAN_UP_FILES - # a custom list of globs, e.g. *.log *.vtk - *.vtk -) -# Usually, you will not need to modify anything beyond this point... - -CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) - -FIND_PACKAGE(deal.II 9.2.0 QUIET - HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} - ) -IF(NOT ${deal.II_FOUND}) - MESSAGE(FATAL_ERROR "\n" - "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n" - "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" - "or set an environment variable \"DEAL_II_DIR\" that contains this path." - ) -ENDIF() - -DEAL_II_INITIALIZE_CACHED_VARIABLES() - -# Store current source directory -SET(TMP ${CMAKE_SOURCE_DIR}) -SET(CMAKE_SOURCE_DIR "..") - -# Query git information -DEAL_II_QUERY_GIT_INFORMATION() - -# Restore source dir -SET(CMAKE_SOURCE_DIR ${TMP}) - -PROJECT(${TARGET} LANGUAGES CXX) - -DEAL_II_INVOKE_AUTOPILOT() - -FIND_PACKAGE(precice REQUIRED) -TARGET_LINK_LIBRARIES(${TARGET} precice::precice) - -SET_PROPERTY(TARGET ${TARGET} APPEND PROPERTY COMPILE_DEFINITIONS - GIT_BRANCH="${GIT_BRANCH}" - GIT_REVISION="${GIT_REVISION}" - GIT_SHORTREV="${GIT_SHORTREV}") diff --git a/nonlinear_elasticity/include/parameter_handling.h b/nonlinear_elasticity/include/parameter_handling.h deleted file mode 100644 index bdb5952..0000000 --- a/nonlinear_elasticity/include/parameter_handling.h +++ /dev/null @@ -1,414 +0,0 @@ -#ifndef parameter_handling_h -#define parameter_handling_h - -#include - -namespace Nonlinear_Elasticity -{ - using namespace dealii; - - /** - * This class declares all parameters, which can be specified in the parameter - * file. The subsection abut preCICE configurations is directly interlinked - * to the Adapter class. - */ - namespace Parameters - { - /** - * @brief System: Specifies system properties - */ - struct System - { - double nu; - double mu; - double rho; - Tensor<1, 3, double> body_force; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - System::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - prm.declare_entry("Poisson's ratio", - "0.3", - Patterns::Double(-1.0, 0.5), - "Poisson's ratio"); - - prm.declare_entry("Shear modulus", - "0.4225e6", - Patterns::Double(), - "Shear modulus"); - - prm.declare_entry("rho", "1000", Patterns::Double(), "rho"); - - prm.declare_entry("body forces", - "0,0,0", - Patterns::List(Patterns::Double()), - "body forces x,y,z"); - } - prm.leave_subsection(); - } - - void - System::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("System properties"); - { - nu = prm.get_double("Poisson's ratio"); - mu = prm.get_double("Shear modulus"); - rho = prm.get_double("rho"); - const std::vector body_forces_input = - Utilities::split_string_list(prm.get("body forces")); - for (uint d = 0; d < 3; ++d) - body_force[d] = Utilities::string_to_double(body_forces_input[d]); - } - prm.leave_subsection(); - } - - /** - * @brief LinearSolver: Specifies linear solver properties - */ - struct LinearSolver - { - std::string type_lin; - double tol_lin; - double max_iterations_lin; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - LinearSolver::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - prm.declare_entry("Solver type", - "CG", - Patterns::Selection("CG|Direct"), - "Linear solver: CG or Direct"); - - prm.declare_entry("Residual", - "1e-6", - Patterns::Double(0.0), - "Linear solver residual (scaled by residual norm)"); - - prm.declare_entry( - "Max iteration multiplier", - "1", - Patterns::Double(0.0), - "Linear solver iterations (multiples of the system matrix size)"); - } - prm.leave_subsection(); - } - - void - LinearSolver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Linear solver"); - { - type_lin = prm.get("Solver type"); - tol_lin = prm.get_double("Residual"); - max_iterations_lin = prm.get_double("Max iteration multiplier"); - } - prm.leave_subsection(); - } - - /** - * @brief NonlinearSolver: Specifies nonlinear solver properties - */ - struct NonlinearSolver - { - unsigned int max_iterations_NR; - double tol_f; - double tol_u; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - NonlinearSolver::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Nonlinear solver"); - { - prm.declare_entry("Max iterations Newton-Raphson", - "10", - Patterns::Integer(0), - "Number of Newton-Raphson iterations allowed"); - - prm.declare_entry("Tolerance force", - "1.0e-9", - Patterns::Double(0.0), - "Force residual tolerance"); - - prm.declare_entry("Tolerance displacement", - "1.0e-6", - Patterns::Double(0.0), - "Displacement error tolerance"); - } - prm.leave_subsection(); - } - - void - NonlinearSolver::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Nonlinear solver"); - { - max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson"); - tol_f = prm.get_double("Tolerance force"); - tol_u = prm.get_double("Tolerance displacement"); - } - prm.leave_subsection(); - } - - /** - * @brief Time: Specifies simulation time properties - */ - struct Time - { - double delta_t; - double end_time; - int output_interval; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - Time::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - prm.declare_entry("End time", "1", Patterns::Double(), "End time"); - - prm.declare_entry("Time step size", - "0.1", - Patterns::Double(), - "Time step size"); - - prm.declare_entry("Output interval", - "1", - Patterns::Integer(), - "Output interval"); - } - prm.leave_subsection(); - } - - void - Time::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Time"); - { - end_time = prm.get_double("End time"); - delta_t = prm.get_double("Time step size"); - output_interval = prm.get_integer("Output interval"); - } - prm.leave_subsection(); - } - - /** - * @brief Discretization: Specifies parameters for time integration by an - * implicit Newmark scheme and polynomial degree of the FE system - */ - struct Discretization - { - double beta; - double gamma; - unsigned int poly_degree; - - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - void - Discretization::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - prm.declare_entry("beta", - "0.25", - Patterns::Double(0, 0.5), - "Newmark beta"); - - prm.declare_entry("gamma", - "0.5", - Patterns::Double(0, 1), - "Newmark gamma"); - - prm.declare_entry("Polynomial degree", - "3", - Patterns::Integer(0), - "Polynomial degree of the FE system"); - } - prm.leave_subsection(); - } - - void - Discretization::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("Discretization"); - { - beta = prm.get_double("beta"); - gamma = prm.get_double("gamma"); - poly_degree = prm.get_integer("Polynomial degree"); - } - prm.leave_subsection(); - } - - /** - * @brief PreciceAdapterConfiguration: Specifies preCICE related information. - * A lot of these information need to be consistent with the - * precice-config.xml file. - */ - struct PreciceAdapterConfiguration - { - std::string scenario; - std::string config_file; - std::string participant_name; - std::string mesh_name; - std::string read_data_name; - std::string write_data_name; - double flap_location; - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - - void - PreciceAdapterConfiguration::declare_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - prm.declare_entry("Scenario", - "FSI3", - Patterns::Selection("FSI3|PF"), - "Cases: FSI3 or PF for perpendicular flap"); - prm.declare_entry("precice config-file", - "precice-config.xml", - Patterns::Anything(), - "Name of the precice configuration file"); - prm.declare_entry( - "Participant name", - "dealiisolver", - Patterns::Anything(), - "Name of the participant in the precice-config.xml file"); - prm.declare_entry( - "Mesh name", - "dealii-mesh", - Patterns::Anything(), - "Name of the coupling mesh in the precice-config.xml file"); - prm.declare_entry( - "Read data name", - "received-data", - Patterns::Anything(), - "Name of the read data in the precice-config.xml file"); - prm.declare_entry( - "Write data name", - "calculated-data", - Patterns::Anything(), - "Name of the write data in the precice-config.xml file"); - prm.declare_entry("Flap location", - "0.0", - Patterns::Double(-3, 3), - "PF x-location"); - } - prm.leave_subsection(); - } - - void - PreciceAdapterConfiguration::parse_parameters(ParameterHandler &prm) - { - prm.enter_subsection("precice configuration"); - { - scenario = prm.get("Scenario"); - config_file = prm.get("precice config-file"); - participant_name = prm.get("Participant name"); - mesh_name = prm.get("Mesh name"); - read_data_name = prm.get("Read data name"); - write_data_name = prm.get("Write data name"); - flap_location = prm.get_double("Flap location"); - } - prm.leave_subsection(); - } - - - - struct AllParameters : public System, - public LinearSolver, - public NonlinearSolver, - public Time, - public Discretization, - public PreciceAdapterConfiguration - - { - AllParameters(const std::string &input_file); - - static void - declare_parameters(ParameterHandler &prm); - - void - parse_parameters(ParameterHandler &prm); - }; - - AllParameters::AllParameters(const std::string &input_file) - { - ParameterHandler prm; - declare_parameters(prm); - prm.parse_input(input_file); - parse_parameters(prm); - - // Optional, if we want to print all parameters in the beginning of the - // simulation - // prm.print_parameters(std::cout,ParameterHandler::Text); - } - - void - AllParameters::declare_parameters(ParameterHandler &prm) - { - System::declare_parameters(prm); - LinearSolver::declare_parameters(prm); - NonlinearSolver::declare_parameters(prm); - Time::declare_parameters(prm); - Discretization::declare_parameters(prm); - PreciceAdapterConfiguration::declare_parameters(prm); - } - - void - AllParameters::parse_parameters(ParameterHandler &prm) - { - System::parse_parameters(prm); - LinearSolver::parse_parameters(prm); - NonlinearSolver::parse_parameters(prm); - Time::parse_parameters(prm); - Discretization::parse_parameters(prm); - PreciceAdapterConfiguration::parse_parameters(prm); - } - } // namespace Parameters -} // namespace Nonlinear_Elasticity - -#endif diff --git a/linear_elasticity/linear_elasticity.prm b/parameters.prm similarity index 53% rename from linear_elasticity/linear_elasticity.prm rename to parameters.prm index 62edcbe..d767d85 100644 --- a/linear_elasticity/linear_elasticity.prm +++ b/parameters.prm @@ -14,37 +14,47 @@ subsection Time end subsection Discretization - # Time integration scheme - # 0 = forward, 1 = backward - set theta = 0.5 - # Polynomial degree of the FE system set Polynomial degree = 3 end subsection System properties - # mu (shear modulus) - set mu = 0.5e6 + # Poisson's ratio + set Poisson's ratio = 0.4 - # lambda - set lambda = 2e6 + # Shear modulus + set Shear modulus = 0.5e6 - # density - set rho = 1000 + # Density + set rho = 1000 - # body forces x,y,z + # Body forces x,y,z set body forces = 0.0,0.0,0.0 end -subsection Linear solver - # Linear solver iterations (multiples of the system matrix size) +subsection Solver + # Structural model to be used: linear or neo-Hookean + set Model = linear + + # Linear solver: CG or Direct + set Solver type = Direct + + # Max CG solver iterations (multiples of the system matrix size) + # In 2D, this value is best set at 2. In 3D, a value of 1 works fine. set Max iteration multiplier = 1 - # Linear solver residual (scaled by residual norm) + # Relative drop criterion for CG solver residual (multiplied by residual norm, ignored if Model == linear or solver == direct) + # Hard-coded to absolute criterion of 1e-10 for the linear model. set Residual = 1e-6 - # Linear solver: CG or Direct - set Solver type = Direct + # Number of Newton-Raphson iterations allowed (ignored if Model == linear) + set Max iterations Newton-Raphson = 10 + + # Relative displacement error tolerance for non-linear iteration (ignored if Model == linear) + set Tolerance displacement = 1.0e-6 + + # Relative force residual tolerance for non-linear iteration (ignored if Model == linear) + set Tolerance force = 1.0e-9 end subsection precice configuration @@ -58,7 +68,7 @@ subsection precice configuration set Participant name = Solid # Name of the coupling mesh in the precice-config.xml file - set Mesh name = Solid_mesh + set Mesh name = Solid-Mesh # Name of the read data in the precice-config.xml file set Read data name = Stress diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt new file mode 100644 index 0000000..d70fb87 --- /dev/null +++ b/source/CMakeLists.txt @@ -0,0 +1,4 @@ +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +ADD_SUBDIRECTORY(nonlinear_elasticity) +ADD_SUBDIRECTORY(linear_elasticity) diff --git a/source/linear_elasticity/CMakeLists.txt b/source/linear_elasticity/CMakeLists.txt new file mode 100644 index 0000000..b5c440b --- /dev/null +++ b/source/linear_elasticity/CMakeLists.txt @@ -0,0 +1,18 @@ +## +# CMake script for the dealii-adapter: +## + +# Set the name of the project and target: +SET(TARGET "linear_elasticity") + +SET(_src + ${TARGET}.cc +) + +CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0) + +ADD_LIBRARY(${TARGET} OBJECT ${_src}) +FIND_PACKAGE(precice REQUIRED) +TARGET_LINK_LIBRARIES(${TARGET} precice::precice) + +DEAL_II_SETUP_TARGET(${TARGET}) diff --git a/source/linear_elasticity/include/linear_elasticity.h b/source/linear_elasticity/include/linear_elasticity.h new file mode 100644 index 0000000..ae270fa --- /dev/null +++ b/source/linear_elasticity/include/linear_elasticity.h @@ -0,0 +1,156 @@ +#ifndef LINEAR_ELASTICITY_H +#define LINEAR_ELASTICITY_H + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "postprocessor.h" + +// The Linear_Elasticity case includes a linear elastic material with a one-step +// theta time integration +namespace Linear_Elasticity +{ + using namespace dealii; + + template + class ElastoDynamics + { + public: + ElastoDynamics(const std::string &case_path, + const std::string ¶meter_file); + + ~ElastoDynamics(); + // As usual in dealii, the run function covers the main time loop of the + // system + void + run(); + + private: + // Create the mesh and set boundary IDs for different boundary conditions + void + make_grid(); + + // Set up the FE system and allocate data structures + void + setup_system(); + + // Compute time invariant matrices e.g. stiffness matrix and mass matrix + void + assemble_system(); + + // Assemble the Neumann contribution i.e. the coupling data obtained from + // the Fluid participant + void + assemble_rhs(); + + void + assemble_consistent_loading(); + + // Solve the linear system + void + solve(); + + // Update the displacement according to the theta scheme + void + update_displacement(); + + // Output results to vtk files + void + output_results() const; + + // Paramter class parsing all user specific input parameters + const Parameters::AllParameters parameters; + + // Boundary IDs, reserved for the respectve application + unsigned int clamped_mesh_id; + unsigned int out_of_plane_clamped_mesh_id; + const unsigned int interface_boundary_id; + + // Dealii typical objects + Triangulation triangulation; + DoFHandler dof_handler; + FESystem fe; + MappingQGeneric mapping; + const unsigned int quad_order; + + AffineConstraints hanging_node_constraints; + + // Matrices used during computations + SparsityPattern sparsity_pattern; + SparseMatrix mass_matrix; + SparseMatrix stiffness_matrix; + SparseMatrix system_matrix; + SparseMatrix stepping_matrix; + + // Time dependent variables + Vector old_velocity; + Vector velocity; + Vector old_displacement; + Vector displacement; + Vector old_stress; + Vector stress; + Vector system_rhs; + + // Body forces e.g. gravity. Values are specified in the input file + const bool body_force_enabled; + Vector body_force_vector; + + // In order to measure some timings + mutable TimerOutput timer; + + // The main adapter objects: The time class keeps track of the current time + // and time steps. The Adapter class includes all functionalities for + // coupling via preCICE. Look at the documentation of the class for more + // information. + Adapter::Time time; + Adapter::Adapter, Parameters::AllParameters> adapter; + + // Alias for all time dependent variables, which should be saved/reloaded + // in case of an implicit coupling. This vector is directly used in the + // Adapter class + std::vector *> state_variables; + // for the output directory + const std::string case_path; + }; +} // namespace Linear_Elasticity + +#endif // LINEAR_ELASTICITY_H diff --git a/linear_elasticity/include/postprocessor.h b/source/linear_elasticity/include/postprocessor.h similarity index 98% rename from linear_elasticity/include/postprocessor.h rename to source/linear_elasticity/include/postprocessor.h index c3abd58..a36b024 100644 --- a/linear_elasticity/include/postprocessor.h +++ b/source/linear_elasticity/include/postprocessor.h @@ -1,5 +1,4 @@ -#ifndef POSTPROCESSOR_H -#define POSTPROCESSOR_H +#pragma once #include @@ -122,4 +121,3 @@ namespace Linear_Elasticity } } // namespace Linear_Elasticity -#endif // POSTPROCESSOR_H diff --git a/linear_elasticity/linear_elasticity.cc b/source/linear_elasticity/linear_elasticity.cc similarity index 81% rename from linear_elasticity/linear_elasticity.cc rename to source/linear_elasticity/linear_elasticity.cc index 6eac8e9..1c7dfc8 100644 --- a/linear_elasticity/linear_elasticity.cc +++ b/source/linear_elasticity/linear_elasticity.cc @@ -1,3 +1,5 @@ +#include "include/linear_elasticity.h" + #include #include #include @@ -33,12 +35,13 @@ #include #include +#include +#include +#include + #include #include -#include "../adapter/adapter.h" -#include "../adapter/time.h" -#include "include/parameter_handling.h" #include "include/postprocessor.h" // The Linear_Elasticity case includes a linear elastic material with a one-step @@ -47,111 +50,11 @@ namespace Linear_Elasticity { using namespace dealii; - template - class ElastoDynamics - { - public: - ElastoDynamics(const std::string &case_path); - ~ElastoDynamics(); - // As usual in dealii, the run function covers the main time loop of the - // system - void - run(); - - private: - // Create the mesh and set boundary IDs for different boundary conditions - void - make_grid(); - - // Set up the FE system and allocate data structures - void - setup_system(); - - // Compute time invariant matrices e.g. stiffness matrix and mass matrix - void - assemble_system(); - - // Assemble the Neumann contribution i.e. the coupling data obtained from - // the Fluid participant - void - assemble_rhs(); - - void - assemble_consistent_loading(); - - // Solve the linear system - void - solve(); - - // Update the displacement according to the theta scheme - void - update_displacement(); - - // Output results to vtk files - void - output_results() const; - - // Paramter class parsing all user specific input parameters - const Parameters::AllParameters parameters; - - // Boundary IDs, reserved for the respectve application - unsigned int clamped_mesh_id; - unsigned int out_of_plane_clamped_mesh_id; - const unsigned int interface_boundary_id; - - // Dealii typical objects - Triangulation triangulation; - DoFHandler dof_handler; - FESystem fe; - MappingQGeneric mapping; - const unsigned int quad_order; - - AffineConstraints hanging_node_constraints; - - // Matrices used during computations - SparsityPattern sparsity_pattern; - SparseMatrix mass_matrix; - SparseMatrix stiffness_matrix; - SparseMatrix system_matrix; - SparseMatrix stepping_matrix; - - // Time dependent variables - Vector old_velocity; - Vector velocity; - Vector old_displacement; - Vector displacement; - Vector old_stress; - Vector stress; - Vector system_rhs; - - // Body forces e.g. gravity. Values are specified in the input file - const bool body_force_enabled; - Vector body_force_vector; - - // In order to measure some timings - mutable TimerOutput timer; - - // The main adapter objects: The time class keeps track of the current time - // and time steps. The Adapter class includes all functionalities for - // coupling via preCICE. Look at the documentation of the class for more - // information. - Adapter::Time time; - Adapter::Adapter, Parameters::AllParameters> adapter; - - // Alias for all time dependent variables, which should be saved/reloaded - // in case of an implicit coupling. This vector is directly used in the - // Adapter class - std::vector *> state_variables; - // for the output directory - const std::string case_path; - }; - - - // Constructor template - ElastoDynamics::ElastoDynamics(const std::string &case_path) - : parameters(case_path + "linear_elasticity.prm") + ElastoDynamics::ElastoDynamics(const std::string &case_path, + const std::string ¶meter_file) + : parameters(parameter_file) , interface_boundary_id(6) , dof_handler(triangulation) , fe(FE_Q(parameters.poly_degree), dim) @@ -638,7 +541,7 @@ namespace Linear_Elasticity const int solver_its = system_matrix.m() * parameters.max_iterations_lin; - const double tol_sol = parameters.tol_lin * system_rhs.l2_norm(); + const double tol_sol = 1.e-10; SolverControl solver_control(solver_its, tol_sol); GrowingVectorMemory<> GVM; @@ -807,82 +710,6 @@ namespace Linear_Elasticity // communication etc. adapter.precice.finalize(); } -} // namespace Linear_Elasticity - -int -main(int argc, char **argv) -{ - using namespace Linear_Elasticity; - using namespace dealii; -#ifdef DEAL_II_WITH_MPI - Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); -#endif - - try - { - // Query adapter and deal.II info - const std::string adapter_info = - GIT_SHORTREV == std::string("") ? - "unknown" : - (GIT_SHORTREV + std::string(" on branch ") + GIT_BRANCH); - const std::string dealii_info = - DEAL_II_GIT_SHORTREV == std::string("") ? - "unknown" : - (DEAL_II_GIT_SHORTREV + std::string(" on branch ") + - DEAL_II_GIT_BRANCH); - - std::cout - << "-----------------------------------------------------------------------------" - << std::endl; - std::cout << "-- . adapter revision " << adapter_info << std::endl; - std::cout << "-- . deal.II " << DEAL_II_PACKAGE_VERSION - << " (revision " << dealii_info << ")" << std::endl; - std::cout - << "-----------------------------------------------------------------------------" - << std::endl - << std::endl; - - std::string parameter_file; - if (argc > 1) - parameter_file = argv[1]; - else - parameter_file = "linear_elasticity.prm"; - - // Extract case path for the output directory - size_t pos = parameter_file.find_last_of("/"); - std::string case_path = - std::string::npos == pos ? "" : parameter_file.substr(0, pos + 1); - - ElastoDynamics elastic_solver(case_path); - elastic_solver.run(); - } - catch (std::exception &exc) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -} + template class ElastoDynamics; +} // namespace Linear_Elasticity diff --git a/source/nonlinear_elasticity/CMakeLists.txt b/source/nonlinear_elasticity/CMakeLists.txt new file mode 100644 index 0000000..3c8a2a5 --- /dev/null +++ b/source/nonlinear_elasticity/CMakeLists.txt @@ -0,0 +1,18 @@ +## +# CMake script for the dealii-adapter: +## + +# Set the name of the project and target: +SET(TARGET "nonlinear_elasticity") + +SET(_src + ${TARGET}.cc +) + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12) + +ADD_LIBRARY(${TARGET} OBJECT ${_src}) +FIND_PACKAGE(precice REQUIRED) +TARGET_LINK_LIBRARIES(${TARGET} precice::precice) + +DEAL_II_SETUP_TARGET(${TARGET}) diff --git a/nonlinear_elasticity/include/compressible_neo_hook_material.h b/source/nonlinear_elasticity/include/compressible_neo_hook_material.h similarity index 97% rename from nonlinear_elasticity/include/compressible_neo_hook_material.h rename to source/nonlinear_elasticity/include/compressible_neo_hook_material.h index ca841aa..75ca8f3 100644 --- a/nonlinear_elasticity/include/compressible_neo_hook_material.h +++ b/source/nonlinear_elasticity/include/compressible_neo_hook_material.h @@ -1,5 +1,4 @@ -#ifndef COMPRESSIBLE_NEO_HOOK_MATERIAL_H -#define COMPRESSIBLE_NEO_HOOK_MATERIAL_H +#pragma once #include @@ -139,4 +138,3 @@ namespace Nonlinear_Elasticity } }; } // namespace Nonlinear_Elasticity -#endif // COMPRESSIBLE_NEO_HOOK_MATERIAL_H diff --git a/source/nonlinear_elasticity/include/nonlinear_elasticity.h b/source/nonlinear_elasticity/include/nonlinear_elasticity.h new file mode 100644 index 0000000..b354c8f --- /dev/null +++ b/source/nonlinear_elasticity/include/nonlinear_elasticity.h @@ -0,0 +1,338 @@ +#ifndef NONLINEAR_ELASTICITY_H +#define NONLINEAR_ELASTICITY_H + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "compressible_neo_hook_material.h" +#include "postprocessor.h" + +namespace Nonlinear_Elasticity +{ + using namespace dealii; + + // PointHistory class offers a method for storing data at the quadrature + // points. Here each quadrature point holds a pointer to a material + // description + template + class PointHistory + { + public: + PointHistory() + {} + + virtual ~PointHistory() + {} + + void + setup_lqp(const Parameters::AllParameters ¶meters) + { + material.reset( + new Material_Compressible_Neo_Hook_One_Field( + parameters.mu, parameters.nu, parameters.rho)); + } + + // Strain energy + NumberType + get_Psi(const NumberType & det_F, + const SymmetricTensor<2, dim, NumberType> &b_bar) const + { + return material->get_Psi(det_F, b_bar); + } + // Kirchhoff stress + SymmetricTensor<2, dim, NumberType> + get_tau(const NumberType & det_F, + const SymmetricTensor<2, dim, NumberType> &b_bar) const + { + return material->get_tau(det_F, b_bar); + } + // Tangent + SymmetricTensor<4, dim, NumberType> + get_Jc(const NumberType & det_F, + const SymmetricTensor<2, dim, NumberType> &b_bar) const + { + return material->get_Jc(det_F, b_bar); + } + // Density + NumberType + get_rho() const + { + return material->get_rho(); + } + + private: + std::shared_ptr> + material; + }; + + + // Forward declarations for classes that will perform assembly of the + // linearized system + template + struct Assembler_Base; + template + struct Assembler; + + // The Solid class is the central class in that it represents the problem at + // hand. It follows the usual scheme in that all it really has is a + // constructor, destructor and a run() function that dispatches all the work + // to private functions of this class: + template + class Solid + { + public: + Solid(const std::string &case_path, const std::string ¶meter_file); + + virtual ~Solid(); + + void + run(); + + private: + // Generates grid and sets boundary IDs, which are needed for different BCs + void + make_grid(); + // Set up the finite element system to be solved + void + system_setup(); + + // Several functions to assemble the system and right hand side matrices + // using multithreading. Each of them comes as a wrapper function, one that + // is executed to do the work in the WorkStream model on one cell, and one + // that copies the work done on this one cell into the global object that + // represents it + void + assemble_system(const BlockVector &solution_delta, + const BlockVector &acceleration, + const BlockVector &external_stress); + + // We use a separate data structure to perform the assembly. It needs access + // to some low-level data, so we simply befriend the class instead of + // creating a complex interface to provide access as necessary. + friend struct Assembler_Base; + friend struct Assembler; + + // Apply Dirichlet boundary conditions on the displacement field + void + make_constraints(const int &it_nr); + + // Create and update the quadrature points. Here, no data needs to be copied + // into a global object, so the copy_local_to_global function is empty: + void + setup_qph(); + + // Solve for the displacement using a Newton-Raphson method. + void + solve_nonlinear_timestep(BlockVector &solution_delta); + + std::pair + solve_linear_system(BlockVector &newton_update); + + // Solution retrieval + BlockVector + get_total_solution(const BlockVector &solution_delta) const; + + // Update fnuctions for time dependent variables according to Newmarks + // scheme + void + update_acceleration(BlockVector &displacement_delta); + + void + update_velocity(BlockVector &displacement_delta); + + void + update_old_variables(); + + // Post-processing and writing data to file : + void + output_results() const; + + const Parameters::AllParameters parameters; + + double vol_reference; + double vol_current; + + Triangulation triangulation; + + CellDataStorage::cell_iterator, + PointHistory> + quadrature_point_history; + + const unsigned int degree; + const FESystem fe; + DoFHandler dof_handler_ref; + const unsigned int dofs_per_cell; + const FEValuesExtractors::Vector u_fe; + + // Description of how the block-system is arranged. There is just 1 block, + // that contains a vector DOF u. This is a legacy of the original work + // (step-44) + static const unsigned int n_blocks = 1; + static const unsigned int n_components = dim; + static const unsigned int first_u_component = 0; + + enum + { + u_dof = 0 + }; + + std::vector dofs_per_block; + + const QGauss qf_cell; + const QGauss qf_face; + const unsigned int n_q_points; + const unsigned int n_q_points_f; + // Interface ID, which is later assigned to the mesh region for coupling + // It is chosen arbotrarily + const unsigned int boundary_interface_id; + + // Newmark parameters + // Coefficients, which are needed for time dependencies + const double alpha_1 = + 1. / (parameters.beta * std::pow(parameters.delta_t, 2)); + const double alpha_2 = 1. / (parameters.beta * parameters.delta_t); + const double alpha_3 = (1 - (2 * parameters.beta)) / (2 * parameters.beta); + const double alpha_4 = + parameters.gamma / (parameters.beta * parameters.delta_t); + const double alpha_5 = 1 - (parameters.gamma / parameters.beta); + const double alpha_6 = + (1 - (parameters.gamma / (2 * parameters.beta))) * parameters.delta_t; + + // Read body force from parameter file + const Tensor<1, 3, double> body_force = parameters.body_force; + + // Clamped boundary ID to be used consistently + const unsigned int clamped_boundary_id = 1; + const unsigned int out_of_plane_clamped_mesh_id = 8; + + // ..and store the directory, in order to output the result files there + const std::string case_path; + + AffineConstraints constraints; + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix tangent_matrix; + BlockVector system_rhs; + BlockVector total_displacement; + BlockVector total_displacement_old; + BlockVector velocity; + BlockVector velocity_old; + BlockVector acceleration; + BlockVector acceleration_old; + + // Alias to collect all time dependent variables in a single vector + // This is directly passed to the Adapter routine in order to + // store these variables for implicit couplings. + std::vector *> state_variables; + + // Global vector, which keeps all contributions of the Fluid participant + // i.e. stress data for assembly. This vector is filled properly in the + // Adapter + BlockVector external_stress; + + // In order to measure some timings + mutable TimerOutput timer; + + // The main adapter objects: The time class keeps track of the current time + // and time steps. The Adapter class includes all functionalities for + // coupling via preCICE. Look at the documentation of the class for more + // information. + Adapter::Time time; + Adapter::Adapter, Parameters::AllParameters> + adapter; + + // Then define a number of variables to store norms and update norms and + // normalisation factors. + struct Errors + { + Errors() + : u(1.0) + {} + + void + reset() + { + u = 1.0; + } + void + normalise(const Errors &val) + { + if (val.u != 0.0) + u /= val.u; + } + + double u; + }; + + Errors error_residual, error_residual_0, error_residual_norm, error_update, + error_update_0, error_update_norm; + + // Methods to calculate erros + void + get_error_residual(Errors &error_residual); + + void + get_error_update(const BlockVector &newton_update, + Errors & error_update); + + // Print information to screen during simulation + static void + print_conv_header(); + + void + print_conv_footer(); + }; + +} // namespace Nonlinear_Elasticity + +#endif // NONLINEAR_ELASTICITY_H diff --git a/nonlinear_elasticity/include/postprocessor.h b/source/nonlinear_elasticity/include/postprocessor.h similarity index 98% rename from nonlinear_elasticity/include/postprocessor.h rename to source/nonlinear_elasticity/include/postprocessor.h index 631b20c..d4acf42 100644 --- a/nonlinear_elasticity/include/postprocessor.h +++ b/source/nonlinear_elasticity/include/postprocessor.h @@ -1,5 +1,4 @@ -#ifndef POSTPROCESSOR_H -#define POSTPROCESSOR_H +#pragma once #include @@ -122,4 +121,3 @@ namespace Nonlinear_Elasticity } } // namespace Nonlinear_Elasticity -#endif // POSTPROCESSOR_H diff --git a/nonlinear_elasticity/nonlinear_elasticity.cc b/source/nonlinear_elasticity/nonlinear_elasticity.cc similarity index 79% rename from nonlinear_elasticity/nonlinear_elasticity.cc rename to source/nonlinear_elasticity/nonlinear_elasticity.cc index 4a20e5c..5821989 100644 --- a/nonlinear_elasticity/nonlinear_elasticity.cc +++ b/source/nonlinear_elasticity/nonlinear_elasticity.cc @@ -1,3 +1,5 @@ +#include "include/nonlinear_elasticity.h" + #include #include #include @@ -45,297 +47,22 @@ #include #include +#include + #include #include -#include "../adapter/adapter.h" -#include "../adapter/time.h" -#include "include/compressible_neo_hook_material.h" -#include "include/parameter_handling.h" #include "include/postprocessor.h" namespace Nonlinear_Elasticity { using namespace dealii; - - // PointHistory class offers a method for storing data at the quadrature - // points. Here each quadrature point holds a pointer to a material - // description - template - class PointHistory - { - public: - PointHistory() - {} - - virtual ~PointHistory() - {} - - void - setup_lqp(const Parameters::AllParameters ¶meters) - { - material.reset( - new Material_Compressible_Neo_Hook_One_Field( - parameters.mu, parameters.nu, parameters.rho)); - } - - // Strain energy - NumberType - get_Psi(const NumberType & det_F, - const SymmetricTensor<2, dim, NumberType> &b_bar) const - { - return material->get_Psi(det_F, b_bar); - } - // Kirchhoff stress - SymmetricTensor<2, dim, NumberType> - get_tau(const NumberType & det_F, - const SymmetricTensor<2, dim, NumberType> &b_bar) const - { - return material->get_tau(det_F, b_bar); - } - // Tangent - SymmetricTensor<4, dim, NumberType> - get_Jc(const NumberType & det_F, - const SymmetricTensor<2, dim, NumberType> &b_bar) const - { - return material->get_Jc(det_F, b_bar); - } - // Density - NumberType - get_rho() const - { - return material->get_rho(); - } - - private: - std::shared_ptr> - material; - }; - - - // Forward declarations for classes that will perform assembly of the - // linearized system - template - struct Assembler_Base; - template - struct Assembler; - - // The Solid class is the central class in that it represents the problem at - // hand. It follows the usual scheme in that all it really has is a - // constructor, destructor and a run() function that dispatches all the work - // to private functions of this class: - template - class Solid - { - public: - Solid(const std::string &case_path); - - virtual ~Solid(); - - void - run(); - - private: - // Generates grid and sets boundary IDs, which are needed for different BCs - void - make_grid(); - // Set up the finite element system to be solved - void - system_setup(); - - // Several functions to assemble the system and right hand side matrices - // using multithreading. Each of them comes as a wrapper function, one that - // is executed to do the work in the WorkStream model on one cell, and one - // that copies the work done on this one cell into the global object that - // represents it - void - assemble_system(const BlockVector &solution_delta, - const BlockVector &acceleration, - const BlockVector &external_stress); - - // We use a separate data structure to perform the assembly. It needs access - // to some low-level data, so we simply befriend the class instead of - // creating a complex interface to provide access as necessary. - friend struct Assembler_Base; - friend struct Assembler; - - // Apply Dirichlet boundary conditions on the displacement field - void - make_constraints(const int &it_nr); - - // Create and update the quadrature points. Here, no data needs to be copied - // into a global object, so the copy_local_to_global function is empty: - void - setup_qph(); - - // Solve for the displacement using a Newton-Raphson method. - void - solve_nonlinear_timestep(BlockVector &solution_delta); - - std::pair - solve_linear_system(BlockVector &newton_update); - - // Solution retrieval - BlockVector - get_total_solution(const BlockVector &solution_delta) const; - - // Update fnuctions for time dependent variables according to Newmarks - // scheme - void - update_acceleration(BlockVector &displacement_delta); - - void - update_velocity(BlockVector &displacement_delta); - - void - update_old_variables(); - - // Post-processing and writing data to file : - void - output_results() const; - - const Parameters::AllParameters parameters; - - double vol_reference; - double vol_current; - - Triangulation triangulation; - - CellDataStorage::cell_iterator, - PointHistory> - quadrature_point_history; - - const unsigned int degree; - const FESystem fe; - DoFHandler dof_handler_ref; - const unsigned int dofs_per_cell; - const FEValuesExtractors::Vector u_fe; - - // Description of how the block-system is arranged. There is just 1 block, - // that contains a vector DOF u. This is a legacy of the original work - // (step-44) - static const unsigned int n_blocks = 1; - static const unsigned int n_components = dim; - static const unsigned int first_u_component = 0; - - enum - { - u_dof = 0 - }; - - std::vector dofs_per_block; - - const QGauss qf_cell; - const QGauss qf_face; - const unsigned int n_q_points; - const unsigned int n_q_points_f; - // Interface ID, which is later assigned to the mesh region for coupling - // It is chosen arbotrarily - const unsigned int boundary_interface_id; - - // Newmark parameters - // Coefficients, which are needed for time dependencies - const double alpha_1 = - 1. / (parameters.beta * std::pow(parameters.delta_t, 2)); - const double alpha_2 = 1. / (parameters.beta * parameters.delta_t); - const double alpha_3 = (1 - (2 * parameters.beta)) / (2 * parameters.beta); - const double alpha_4 = - parameters.gamma / (parameters.beta * parameters.delta_t); - const double alpha_5 = 1 - (parameters.gamma / parameters.beta); - const double alpha_6 = - (1 - (parameters.gamma / (2 * parameters.beta))) * parameters.delta_t; - - // Read body force from parameter file - const Tensor<1, 3, double> body_force = parameters.body_force; - - // Clamped boundary ID to be used consistently - const unsigned int clamped_boundary_id = 1; - const unsigned int out_of_plane_clamped_mesh_id = 8; - - // ..and store the directory, in order to output the result files there - const std::string case_path; - - AffineConstraints constraints; - BlockSparsityPattern sparsity_pattern; - BlockSparseMatrix tangent_matrix; - BlockVector system_rhs; - BlockVector total_displacement; - BlockVector total_displacement_old; - BlockVector velocity; - BlockVector velocity_old; - BlockVector acceleration; - BlockVector acceleration_old; - - // Alias to collect all time dependent variables in a single vector - // This is directly passed to the Adapter routine in order to - // store these variables for implicit couplings. - std::vector *> state_variables; - - // Global vector, which keeps all contributions of the Fluid participant - // i.e. stress data for assembly. This vector is filled properly in the - // Adapter - BlockVector external_stress; - - // In order to measure some timings - mutable TimerOutput timer; - - // The main adapter objects: The time class keeps track of the current time - // and time steps. The Adapter class includes all functionalities for - // coupling via preCICE. Look at the documentation of the class for more - // information. - Adapter::Time time; - Adapter::Adapter, Parameters::AllParameters> - adapter; - - // Then define a number of variables to store norms and update norms and - // normalisation factors. - struct Errors - { - Errors() - : u(1.0) - {} - - void - reset() - { - u = 1.0; - } - void - normalise(const Errors &val) - { - if (val.u != 0.0) - u /= val.u; - } - - double u; - }; - - Errors error_residual, error_residual_0, error_residual_norm, error_update, - error_update_0, error_update_norm; - - // Methods to calculate erros - void - get_error_residual(Errors &error_residual); - - void - get_error_update(const BlockVector &newton_update, - Errors & error_update); - - // Print information to screen during simulation - static void - print_conv_header(); - - void - print_conv_footer(); - }; - - // Constructor initializes member variables and reads the parameter file template - Solid::Solid(const std::string &case_path) - : parameters( - Parameters::AllParameters(case_path + "nonlinear_elasticity.prm")) + Solid::Solid(const std::string &case_path, + const std::string ¶meter_file) + : parameters(Parameters::AllParameters(parameter_file)) , vol_reference(0.0) , vol_current(0.0) , triangulation(Triangulation::maximum_smoothing) @@ -355,7 +82,13 @@ namespace Nonlinear_Elasticity , timer(std::cout, TimerOutput::summary, TimerOutput::wall_times) , time(parameters.end_time, parameters.delta_t) , adapter(parameters, boundary_interface_id) - {} + { + AssertThrow( + parameters.data_consistent, + ExcMessage( + "The neo-Hookean solid doesn't support 'Force' data reading. Please switch to 'Stress' " + "data on the Fluid side or use the linear model of the solid solver")); + } // Destructor clears the DoFHandler template @@ -1499,90 +1232,5 @@ namespace Nonlinear_Elasticity timer.leave_subsection("Output results"); } + template class Solid; } // namespace Nonlinear_Elasticity - -int -main(int argc, char **argv) -{ - using namespace Nonlinear_Elasticity; - using namespace dealii; - -#ifdef DEAL_II_WITH_MPI - Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); -#endif - - try - { - deallog.depth_console(0); - static const unsigned int n_threads = MultithreadInfo::n_threads(); - - // Query adapter and deal.II info - const std::string adapter_info = - GIT_SHORTREV == std::string("") ? - "unknown" : - (GIT_SHORTREV + std::string(" on branch ") + GIT_BRANCH); - const std::string dealii_info = - DEAL_II_GIT_SHORTREV == std::string("") ? - "unknown" : - (DEAL_II_GIT_SHORTREV + std::string(" on branch ") + - DEAL_II_GIT_BRANCH); - - std::cout - << "-----------------------------------------------------------------------------" - << std::endl - << "-- . running with " << n_threads << " thread" - << (n_threads == 1 ? "" : "s") << std::endl; - - std::cout << "-- . adapter revision " << adapter_info << std::endl; - std::cout << "-- . deal.II " << DEAL_II_PACKAGE_VERSION - << " (revision " << dealii_info << ")" << std::endl; - std::cout - << "-----------------------------------------------------------------------------" - << std::endl - << std::endl; - - - std::string parameter_file; - if (argc > 1) - parameter_file = argv[1]; - else - parameter_file = "nonlinear_elasticity.prm"; - - // Extract case path for the output directory - size_t pos = parameter_file.find_last_of("/"); - std::string case_path = - std::string::npos == pos ? "" : parameter_file.substr(0, pos + 1); - - // Dimension is determinded via cmake -DDIM - Solid solid(case_path); - solid.run(); - } - catch (std::exception &exc) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Exception on processing: " << std::endl - << exc.what() << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - - return 1; - } - catch (...) - { - std::cerr << std::endl - << std::endl - << "----------------------------------------------------" - << std::endl; - std::cerr << "Unknown exception!" << std::endl - << "Aborting!" << std::endl - << "----------------------------------------------------" - << std::endl; - return 1; - } - - return 0; -} diff --git a/nonlinear_elasticity/nonlinear_elasticity.prm b/source/nonlinear_elasticity/nonlinear_elasticity.prm similarity index 100% rename from nonlinear_elasticity/nonlinear_elasticity.prm rename to source/nonlinear_elasticity/nonlinear_elasticity.prm