diff --git a/CMakeLists.txt b/CMakeLists.txt index 107e598..43bb5fc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -29,47 +29,22 @@ endif() add_definitions(-DSTRAN_DATA_DIR="${CMAKE_CURRENT_SOURCE_DIR}/filedata/") # add example files -add_custom_target(${PROJECT_NAME}-all) -add_custom_target(${PROJECT_NAME}-examples) -aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR}/examples FILES) -foreach(file ${FILES}) - add_gismo_executable(${file}) - get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") - set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") - add_dependencies(${PROJECT_NAME}-examples ${tarname}) - add_dependencies(${PROJECT_NAME}-all ${tarname}) - # install the example executables (optionally) - install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) -endforeach(file ${FILES}) - -# add benchmark files -add_custom_target(${PROJECT_NAME}-benchmarks) -aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR}/benchmarks FILES) -foreach(file ${FILES}) - add_gismo_executable(${file}) - get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") - set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") - add_dependencies(${PROJECT_NAME}-benchmarks ${tarname}) - add_dependencies(${PROJECT_NAME}-all ${tarname}) - # install the example executables (optionally) - install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) -endforeach(file ${FILES}) - -# add solver files -add_custom_target(${PROJECT_NAME}-solvers) -aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR}/solvers FILES) -foreach(file ${FILES}) - add_gismo_executable(${file}) - get_filename_component(tarname ${file} NAME_WE) # name without extension - set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") - set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") - add_dependencies(${PROJECT_NAME}-solvers ${tarname}) - add_dependencies(${PROJECT_NAME}-all ${tarname}) - # install the example executables (optionally) - install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) -endforeach(file ${FILES}) +if(GISMO_BUILD_EXAMPLES) + add_custom_target(${PROJECT_NAME}-all) + add_custom_target(${PROJECT_NAME}-examples) + add_custom_target(${PROJECT_NAME}-benchmarks) + add_custom_target(${PROJECT_NAME}-solvers) + add_custom_target(${PROJECT_NAME}-tutorials) + add_subdirectory(examples) + add_subdirectory(benchmarks) + add_subdirectory(solvers) + add_subdirectory(tutorials) +else() + add_subdirectory(examples EXCLUDE_FROM_ALL) + add_subdirectory(benchmarks EXCLUDE_FROM_ALL) + add_subdirectory(solvers EXCLUDE_FROM_ALL) + add_subdirectory(tutorials EXCLUDE_FROM_ALL) +endif(GISMO_BUILD_EXAMPLES) # # add unittests # aux_gs_cpp_directory(${PROJECT_SOURCE_DIR}/unittests unittests_SRCS) diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt new file mode 100644 index 0000000..06dcbf9 --- /dev/null +++ b/benchmarks/CMakeLists.txt @@ -0,0 +1,26 @@ +###################################################################### +## CMakeLists.txt --- gsStructuralAnalysis/benchmarks +## This file is part of the G+Smo library. +## +## Author: Angelos Mantzaflaris, Hugo Verhelst +## Copyright (C) 2023 +###################################################################### + +# add benchmark files +aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) +foreach(file ${FILES}) + add_gismo_executable(${file}) + get_filename_component(tarname ${file} NAME_WE) # name without extension + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + if(GISMO_BUILD_EXAMPLES) + set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + else(GISMO_BUILD_EXAMPLES) + set_target_properties(${tarname} PROPERTIES + FOLDER "${PROJECT_NAME}" + EXCLUDE_FROM_ALL TRUE) + endif(GISMO_BUILD_EXAMPLES) + add_dependencies(${PROJECT_NAME}-benchmarks ${tarname}) + add_dependencies(${PROJECT_NAME}-all ${tarname}) + # install the example executables (optionally) + install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) +endforeach(file ${FILES}) \ No newline at end of file diff --git a/benchmarks/benchmark_Balloon.cpp b/benchmarks/benchmark_Balloon.cpp index 2dda304..3ab8f85 100644 --- a/benchmarks/benchmark_Balloon.cpp +++ b/benchmarks/benchmark_Balloon.cpp @@ -152,11 +152,7 @@ int main (int argc, char** argv) dirname = dirname + "/Balloon"; std::string output = "solution"; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); - + gsFileManager::mkdir(dirname); // plot geometry if (plot) diff --git a/benchmarks/benchmark_Beam.cpp b/benchmarks/benchmark_Beam.cpp index 6d393aa..bd4bac1 100644 --- a/benchmarks/benchmark_Beam.cpp +++ b/benchmarks/benchmark_Beam.cpp @@ -213,10 +213,7 @@ int main (int argc, char** argv) else GISMO_ERROR("Testcase unknown"); - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir(dirname); // plot geometry if (plot) diff --git a/benchmarks/benchmark_Cylinder.cpp b/benchmarks/benchmark_Cylinder.cpp index 88f9457..e0fa239 100644 --- a/benchmarks/benchmark_Cylinder.cpp +++ b/benchmarks/benchmark_Cylinder.cpp @@ -148,11 +148,7 @@ int main (int argc, char** argv) std::string output = "solution"; wn = output + "data.txt"; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); - + gsFileManager::mkdir(dirname); // plot geometry if (plot) diff --git a/benchmarks/benchmark_Elasticity_Beam_APALM.cpp b/benchmarks/benchmark_Elasticity_Beam_APALM.cpp index 9a85f6a..95b1934 100644 --- a/benchmarks/benchmark_Elasticity_Beam_APALM.cpp +++ b/benchmarks/benchmark_Elasticity_Beam_APALM.cpp @@ -283,9 +283,7 @@ int main (int argc, char** argv) wn = output + "data.txt"; } - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - system(command); + gsFileManager::mkdir(dirname); gsFunctionExpr<> g("0","0","0",3); diff --git a/benchmarks/benchmark_FrustrumALM.cpp b/benchmarks/benchmark_FrustrumALM.cpp index b55c36e..7976b76 100644 --- a/benchmarks/benchmark_FrustrumALM.cpp +++ b/benchmarks/benchmark_FrustrumALM.cpp @@ -186,10 +186,7 @@ int main (int argc, char** argv) wn = output + "data.txt"; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir(dirname); // plot geometry if (plot) diff --git a/benchmarks/benchmark_FrustrumDC.cpp b/benchmarks/benchmark_FrustrumDC.cpp index fcac7df..7ee78f3 100644 --- a/benchmarks/benchmark_FrustrumDC.cpp +++ b/benchmarks/benchmark_FrustrumDC.cpp @@ -176,10 +176,7 @@ int main (int argc, char** argv) wn = output + "data.txt"; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir(dirname); // plot geometry if (plot) diff --git a/benchmarks/benchmark_MaterialTest.cpp b/benchmarks/benchmark_MaterialTest.cpp index 8717861..55a8e9b 100644 --- a/benchmarks/benchmark_MaterialTest.cpp +++ b/benchmarks/benchmark_MaterialTest.cpp @@ -277,10 +277,7 @@ int main (int argc, char** argv) output = "solution"; wn = output + "data.txt"; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir(dirname); // plot geometry if (plot) diff --git a/benchmarks/benchmark_MaterialTestConv.cpp b/benchmarks/benchmark_MaterialTestConv.cpp index 39323e8..a3ec660 100644 --- a/benchmarks/benchmark_MaterialTestConv.cpp +++ b/benchmarks/benchmark_MaterialTestConv.cpp @@ -306,7 +306,7 @@ int main (int argc, char** argv) staticSolver.setOptions(solverOptions); gsStatus status = staticSolver.solveNonlinear(solVector); - GISMO_ASSERT(status==gsStatus::Success,"Newton solver failed"); + GISMO_ENSURE(status==gsStatus::Success,"Newton solver failed"); mp_def = assembler->constructSolution(solVector); gsMultiPatch<> deformation = mp_def; diff --git a/benchmarks/benchmark_Roof.cpp b/benchmarks/benchmark_Roof.cpp index 7fcd9e1..df7e0b5 100644 --- a/benchmarks/benchmark_Roof.cpp +++ b/benchmarks/benchmark_Roof.cpp @@ -241,10 +241,7 @@ int main (int argc, char** argv) output = "solution"; wn = output + "data.txt"; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir(dirname); // plot geometry diff --git a/benchmarks/benchmark_Roof_DWR.cpp b/benchmarks/benchmark_Roof_DWR.cpp index 11d0174..1bf0655 100644 --- a/benchmarks/benchmark_Roof_DWR.cpp +++ b/benchmarks/benchmark_Roof_DWR.cpp @@ -317,10 +317,7 @@ int main (int argc, char** argv) output = "solution"; wn = output + "data.txt"; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir(dirname); // plot geometry diff --git a/benchmarks/benchmark_TensionWrinkling.cpp b/benchmarks/benchmark_TensionWrinkling.cpp index 6a76607..73806e1 100644 --- a/benchmarks/benchmark_TensionWrinkling.cpp +++ b/benchmarks/benchmark_TensionWrinkling.cpp @@ -251,10 +251,7 @@ int main (int argc, char** argv) real_t cross_val = 0.0; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir(dirname); // plot geometry diff --git a/benchmarks/benchmark_UniaxialTension.cpp b/benchmarks/benchmark_UniaxialTension.cpp index 6c7caa9..be3bcad 100644 --- a/benchmarks/benchmark_UniaxialTension.cpp +++ b/benchmarks/benchmark_UniaxialTension.cpp @@ -141,10 +141,7 @@ int main (int argc, char** argv) std::string output = dirname + "/UniaxialTension"; output = "solution"; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir(dirname); // plot geometry diff --git a/benchmarks/benchmark_Wrinkling.cpp b/benchmarks/benchmark_Wrinkling.cpp index c32dffe..bb717e5 100644 --- a/benchmarks/benchmark_Wrinkling.cpp +++ b/benchmarks/benchmark_Wrinkling.cpp @@ -30,9 +30,11 @@ #include #include +#ifdef gsUnstructuredSplines_ENABLED #include #include #include +#endif using namespace gismo; @@ -42,9 +44,10 @@ void initStepOutput( const std::string name, const gsMatrix & points); template void writeStepOutput(const gsALMBase * arcLength, const gsMultiPatch & deformation, const std::string name, const gsMatrix & points, const gsMatrix & patches); -#ifdef gsKLShell_ENABLED int main (int argc, char** argv) { +#ifdef gsKLShell_ENABLED +#ifdef gsUnstructuredSplines_ENABLED // Input options int numElevate = 2; int numHref = 5; @@ -1073,14 +1076,16 @@ int main (int argc, char** argv) delete arcLength; return result; -} + +#else//gsUnstructuredSplines_ENABLED + gsWarn<<"G+Smo is not compiled with the gsUnstructuredSplines module."; + return EXIT_FAILURE; +#endif #else//gsKLShell_ENABLED -int main(int argc, char *argv[]) -{ gsWarn<<"G+Smo is not compiled with the gsKLShell module."; return EXIT_FAILURE; -} #endif +} template void initStepOutput(const std::string name, const gsMatrix & points) diff --git a/benchmarks/benchmark_Wrinkling_DWR.cpp b/benchmarks/benchmark_Wrinkling_DWR.cpp index a81f60c..dff67c8 100644 --- a/benchmarks/benchmark_Wrinkling_DWR.cpp +++ b/benchmarks/benchmark_Wrinkling_DWR.cpp @@ -307,10 +307,7 @@ int main (int argc, char** argv) index_t cross_coordinate = 0; real_t cross_val = 0.0; - - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - system(command); + gsFileManager::mkdir(dirname); // plot geometry if (plot) diff --git a/doc/gsALMBase.dox b/doc/gsALMBase.dox deleted file mode 100644 index d98f8ba..0000000 --- a/doc/gsALMBase.dox +++ /dev/null @@ -1,32 +0,0 @@ -namespace gismo -{ - -/** \page gsALMBase gsALM-classes -\defgroup gsALMBase gsALMBase -\ingroup gsStructuralAnalysis - -Performs the arc length method to solve a nonlinear equation system. - -More details can be found in [1]. - -Examples with the use of arc-length methods are -- \subpage benchmark_Balloon.cpp -- \subpage benchmark_Beam.cpp -- \subpage benchmark_Cylinder.cpp -- \subpage benchmark_FrustrumALM.cpp -- \subpage benchmark_MaterialTest.cpp -- \subpage benchmark_Roof.cpp -- \subpage benchmark_UniaxialTension.cpp -- \subpage benchmark_Wrinkling.cpp -- \subpage gsElasticity_ArcLength.cpp -- \subpage gsThinShell_ArcLength.cpp -- \subpage gsThinShell_Wrinkling.cpp -- \subpage gsThinShell_WrinklingFitting.cpp -- \subpage gsThinShell_WrinklingPerturbed.cpp - - -[1] de Borst, R., Crisfield, M. A., Remmers, J. J. C., & Verhoosel, C. V. (2012). Non-Linear Finite Element Analysis of Solids and Structures: Second Edition. In Non-Linear Finite Element Analysis of Solids and Structures: Second Edition. https://doi.org/10.1002/9781118375938 - -*/ -} - diff --git a/doc/gsALMSolvers.dox b/doc/gsALMSolvers.dox new file mode 100644 index 0000000..54271c8 --- /dev/null +++ b/doc/gsALMSolvers.dox @@ -0,0 +1,31 @@ +namespace gismo +{ + +/** \defgroup gsALMSolvers Arc-length methods +\ingroup gsStructuralAnalysis + +Performs the arc length method to solve a nonlinear equation system. + +More details can be found in [1]. + +Examples with the use of arc-length methods are +- \ref benchmark_Balloon.cpp +- \ref benchmark_Beam.cpp +- \ref benchmark_Cylinder.cpp +- \ref benchmark_FrustrumALM.cpp +- \ref benchmark_MaterialTest.cpp +- \ref benchmark_Roof.cpp +- \ref benchmark_UniaxialTension.cpp +- \ref benchmark_Wrinkling.cpp +- \ref gsElasticity_ArcLength.cpp +- \ref gsThinShell_ArcLength.cpp +- \ref gsThinShell_Wrinkling.cpp +- \ref gsThinShell_WrinklingFitting.cpp +- \ref gsThinShell_WrinklingPerturbed.cpp + + +[1] de Borst, R., Crisfield, M. A., Remmers, J. J. C., & Verhoosel, C. V. (2012). Non-Linear Finite Element Analysis of Solids and Structures: Second Edition. In Non-Linear Finite Element Analysis of Solids and Structures: Second Edition. https://doi.org/10.1002/9781118375938 + +*/ +} + diff --git a/doc/gsBucklingSolver.dox b/doc/gsDynamicSolvers.dox similarity index 84% rename from doc/gsBucklingSolver.dox rename to doc/gsDynamicSolvers.dox index 2dfecb3..26ebb8b 100644 --- a/doc/gsBucklingSolver.dox +++ b/doc/gsDynamicSolvers.dox @@ -1,8 +1,7 @@ namespace gismo { -/** \page gsBucklingSolver gsBucklingSolver -\defgroup gsBucklingSolver gsBucklingSolver +/** \defgroup gsDynamicSolvers Dynamic solvers \ingroup gsStructuralAnalysis Linear buckling analysis is performed by solving the following eigenvalue problem: @@ -18,7 +17,7 @@ obtained by solving a linear problem \f$K_L\mathbf{u}^L_h = \mathbf{P}\f$. Furth The modeshape is represented by \f$\phi\f$. Examples with the use of this class are: -- \subpage gsThinShell_Buckling.cpp +- \ref gsThinShell_Buckling.cpp */ diff --git a/doc/gsEigenSolvers.dox b/doc/gsEigenSolvers.dox new file mode 100644 index 0000000..f1a15a1 --- /dev/null +++ b/doc/gsEigenSolvers.dox @@ -0,0 +1,41 @@ +namespace gismo +{ + +/** \defgroup gsEigenSolvers Eigenproblem solvers +\ingroup gsStructuralAnalysis + +\section gsEigenSolvers_buckling Linear buckling analysis + +Linear buckling analysis is performed by solving the following eigenvalue problem: + +\f[ + (K_L-\sigma K_{NL}(\mathbf{u}^L_h))\mathbf{v}_h = \lambda K_{NL}(\mathbf{u}_h) \mathbf{v}_h +\f] + +Where \f$K_L\f$ is the linear stiffness matrix, \f$K_{NL}(\mathbf{u}_h)\f$ is the tangential +stiffness matrix assembled around \f$\mathbf{u}^L_h\f$. The solution \f$\mathbf{u}^L_h\f$ is +obtained by solving a linear problem \f$K_L\mathbf{u}^L_h = \mathbf{P}\f$. Furthermore, +\f$\sigma\f$ is a shift and \f$(\lambda+\sigma\f)\mathbf{P}\f$ is the critical buckling load. +The modeshape is represented by \f$\phi\f$. + +Examples with the use of this class are: +- \ref gsThinShell_Buckling.cpp + +\section gsEigenSolvers_modal Linear modal analysis + +Linear modal analysis is performed by solving the following eigenvalue problem: + +\f{align*}{ + (K - \sigma M)\mathbf{v}_h = \lambda M\mathbf{v}_h +\f} + +Where \f$K\f$ is the linear stiffness matrix, \f$M\f$ is the mass matrix. Furthermore, +\f$\sigma\f$ is a shift and \f$(\lambda+\sigma\f)\f$ is the eigenfrequency. +The modeshape is represented by \f$\phi\f$. + +Examples with the use of this class are: +- \ref gsThinShell_Modal.cpp + +*/ + +} diff --git a/doc/gsModalSolver.dox b/doc/gsModalSolver.dox deleted file mode 100644 index 54dde67..0000000 --- a/doc/gsModalSolver.dox +++ /dev/null @@ -1,20 +0,0 @@ -namespace gismo -{ - -/** \page gsModalSolver gsModalSolver -\defgroup gsModalSolver gsModalSolver -\ingroup gsStructuralAnalysis - -Linear modal analysis is performed by solving the following eigenvalue problem: - -\f{align*}{ - (K - \sigma M)\mathbf{v}_h = \lambda M\mathbf{v}_h -\f} - -Where \f$K\f$ is the linear stiffness matrix, \f$M\f$ is the mass matrix. Furthermore, -\f$\sigma\f$ is a shift and \f$(\lambda+\sigma\f)\f$ is the eigenfrequency. -The modeshape is represented by \f$\phi\f$. - -Examples with the use of this class are: -- \subpage gsThinShell_Modal.cpp -} diff --git a/doc/gsStaticBase.dox b/doc/gsStaticBase.dox deleted file mode 100644 index 8d256ef..0000000 --- a/doc/gsStaticBase.dox +++ /dev/null @@ -1,19 +0,0 @@ -namespace gismo -{ - -/** \defgroup gsStaticBase gsStaticBase -\ingroup gsStructuralAnalysis - -Base class for static solvers. - -Examples with the use of static solvers are: -- \subpage benchmark_FrustrumDC.cpp -- \subpage benchmark_MaterialTestConv.cpp -- \subpage example_StaticSolvers.cpp -- \subpage gsElasticity_Static.cpp -- \subpage gsElasticity_Static_XML.cpp -- \subpage gsThinShell_Static.cpp -- \subpage gsThinShell_Static_XML.cpp -- \subpage gsThinShell_WrinklingPerturbedDC.cpp -- \subpage gsThinShell_WrinklingPerturbedDR.cpp -} diff --git a/doc/gsStaticSolvers.dox b/doc/gsStaticSolvers.dox new file mode 100644 index 0000000..f5b9b6a --- /dev/null +++ b/doc/gsStaticSolvers.dox @@ -0,0 +1,22 @@ +namespace gismo +{ + +/** \defgroup gsStaticSolvers Static solvers +\ingroup gsStructuralAnalysis + +Base class for static solvers. + +Examples with the use of static solvers are: +- \ref benchmark_FrustrumDC.cpp +- \ref benchmark_MaterialTestConv.cpp +- \ref example_StaticSolvers.cpp +- \ref gsElasticity_Static.cpp +- \ref gsElasticity_Static_XML.cpp +- \ref gsThinShell_Static.cpp +- \ref gsThinShell_Static_XML.cpp +- \ref gsThinShell_WrinklingPerturbedDC.cpp +- \ref gsThinShell_WrinklingPerturbedDR.cpp + +*/ + +} diff --git a/doc/gsStructuralAnalysis.dox b/doc/gsStructuralAnalysis.dox index e7aab67..b02dbbe 100644 --- a/doc/gsStructuralAnalysis.dox +++ b/doc/gsStructuralAnalysis.dox @@ -9,31 +9,30 @@ namespace gismo { -/** \page gsStructuralAnalysis Structural Analysis module -\defgroup gsStructuralAnalysis gsStructuralAnalysis +/** \defgroup gsStructuralAnalysis Structural Analysis Module \ingroup Modules This module contains algorithms for structural analysis. The module is versatile in its use with the elasticity elements from \ref Elasticity and \ref KLShell. The module is enhanced by the \ref gsSpectra class for eigenvalue computations. -\section IntroStructuralAnalysis Using the Structural Analysis module +\section structuralAnalysis_Intro Using the Structural Analysis module -The routines in the structural analysis module are typically based on four entities: +The routines in the structural analysis module are typically based on different operators: - The external force vector (\f$\mathbf{F}_{\text{ext}}\f$), - The linear stiffness matrix (\f$K_L\f$), - The residual vector (\f$\mathbf{R}(\mathbf{u}_h)\f$), - The tangential stiffness matrix or Jacobian (\f$K(\mathbf{u}_h)\f$), -Where \f$\mathbf{u}_h\f$ is the discrete solution. +Where \f$\mathbf{u}_h\f$ is the discrete solution. All types for the operators used in the \ref gsStructuralAnalysis module are defined by \ref gsStructuralAnalysisOps -\subsection forceVector The force vector +\subsection structuralAnalysis_forceVector The force vector The force vector (\f$\mathbf{F}_{\text{ext}}\f$) is the vector that contains the external force contributions on the object. These could be due to body forces, point loads, etc. It is defined as a \ref gsVector throughout the module. -\subsection linearMatrix The linear stiffness matrix +\subsection structuralAnalysis_linearMatrix The linear stiffness matrix The linear stiffness matrix (\f$K_L\f$) is the matrix that has to be assembled given the undeformed configuration of the object. In fact, it is similar to assembling the tangential stiffness matrix around the zero vector (\f$K(\mathbf{0}_h)\f$). The linear stiffness matrix should always be provided as a \ref gsSparseMatrix. -\subsection Residual The tangential stiffness matrix or Jacobian +\subsection structuralAnalysis_Residual The tangential stiffness matrix or Jacobian The residual vector is the vector that represents the difference between the internal and external forces inside the body: \f{align*}{ @@ -73,7 +72,7 @@ ALResidual_t ALResidual = [&time,&stopwatch,&assembler,&mp_def](gsVector }; ~~~~~ -\subsection Jacobian The tangential stiffness matrix or Jacobian +\subsection structuralAnalysis_Jacobian The tangential stiffness matrix or Jacobian The tangential stiffness matrix is the matrix that follows from linearizing the linear stiffness matrix with respect to the deformations collected a discrete vector (\f$\mathbf{u}_h\f$). Therefore, the tangential stiffness matrix is an object that takes the solution vector as an input. In the module, the tangential stiffness matrix has the following type: ~~~~~ @@ -90,6 +89,13 @@ Jacobian_t Jacobian = [&](gsVector }; ~~~~~ + +\subsection structuralAnalysis_Tutorials Tutorials +In the folder \c gsStructuralAnalysis/tutorials, some tutorials are provided, explaining the use of the \ref gsStructuralAnalysis module. These tutorials follow-up on the tutorials provided in the \ref gsKLShell module, see \ref klShell_Tutorials. The tutorials can be compiled using the target \c gsStructuralAnalysis-tutorials. The available tutorials are: +1. \ref nonlinear_shell_static +2. \ref nonlinear_shell_arclength +3. \ref nonlinear_shell_dynamic + */ } diff --git a/doc/gsThinShell_Buckling.dox b/doc/gsThinShell_Buckling.dox index 8b518a4..f02540a 100644 --- a/doc/gsThinShell_Buckling.dox +++ b/doc/gsThinShell_Buckling.dox @@ -686,8 +686,7 @@ int main (int argc, char** argv) if (plot) { gsInfo<<"Plotting in Paraview...\n"; - int systemRet = system("mkdir -p BucklingResults"); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir("BucklingResults"); gsMultiPatch<> deformation = solution; gsMatrix<> modeShape; @@ -730,8 +729,7 @@ int main (int argc, char** argv) if (write) { - int systemRet = system("mkdir -p BucklingResults"); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir("BucklingResults"); std::string wnM = "BucklingResults/eigenvalues.txt"; writeToCSVfile(wnM,values); } diff --git a/doc/gsTimeIntegrator.dox b/doc/gsTimeIntegrator.dox deleted file mode 100644 index 67de03a..0000000 --- a/doc/gsTimeIntegrator.dox +++ /dev/null @@ -1,13 +0,0 @@ -namespace gismo -{ - -/** \defgroup gsTimeIntegrator gsTimeIntegrator -\ingroup gsStructuralAnalysis - -Provides temporal solvers for structural analysis problems - -Examples with the use of this class are: -- \subpage example_DynamicBeam.cpp -- \subpage example_DynamicShell.cpp -- \subpage example_DynamicShellNL.cpp -} diff --git a/doc/nonlinear_shell_arclength.dox b/doc/nonlinear_shell_arclength.dox new file mode 100644 index 0000000..3c27437 --- /dev/null +++ b/doc/nonlinear_shell_arclength.dox @@ -0,0 +1,44 @@ +namespace gismo { +/** + +\page nonlinear_shell_arclength Tutorial: Non-Linear Quasi-Static analysis using Kirchhoff-Love shells + +This example shows how to perform a non-linear quasi-static analysis using the \ref gsALMSolvers submodule of \ref gsStructuralAnalysis. The arc-length methods provided by this submodule work on user-defined operators from \ref gsStructuralAnalysisOps, hence can be used with any element and in any application. In this example, we use the \ref gsKLShell. The goals of the tutorial are the following: +- Solve a non-linear quasi-static problem using \ref gsALMCrisfield + +For detailed information about the construction of \ref gsThinShellAssembler and \ref gsMaterialMatrixBase, please consult the tutorials \ref linear_shell and \ref nonlinear_shell. + +\subsection nonlinear_shell_arclength_Setup Problem set-up +The set-up of the problem to be solved for this example is identical to the one in \ref nonlinear_shell. The reader is referred to that example for a full overview of the code. + +\subsection nonlinear_shell_arclength_Jacobian Jacobian and Residual +To use the Crisfield arc-length method, or any other arc-length method in \ref gsALMSolvers, we need to define a Jacobian matrix and an arc-length residual. Th former is defined from the \ref gsStructuralAnalysisOps, similar to \ref nonlinear_shell_static. The latter is defined as follows: + +\snippet nonlinear_shell_arclength.cpp Define nonlinear residual functions + +Here, the external load vector \f$F\f$ is subtracted with a factor \f$\lambda\f$ to incorporate the load magnification factor in the residual. + +\subsection nonlinear_shell_arclength_Solver Define the arc-length solver +Using the external force vector \f$F\f$, a Jacobian matrix \f$K_{NL}(\mathbf{u})\f$ and a residual vector \f$R(\mathbf{u})\f$, an arc-length solver can be defined. In the following, a \ref gsALMCrisfield is defined: + +\snippet nonlinear_shell_arclength.cpp Set ALM solver + +\subsection nonlinear_shell_arclength_Solve Solve the non-linear quasi-static problem +The quasi-static problem is solved for combinations of the load factor \f$\lambda\f$ and the solution \f$\mathbf{u}\f$, also known as continuation. The full solution loop is given by: + +\snippet nonlinear_shell_arclength.cpp Solve nonlinear problem + +It can be seen that the export of the solution field is similar as in \ref linear_shell \ref nonlinear_shell, but the Paraview files are added to a \ref gsParaviewCollection. This is saved at the end of the simulation as follows: + +\snippet nonlinear_shell_arclength.cpp Save the paraview collection + +\section nonlinear_shell_arclength_Annotated Annotated source file + +Here is the full file \c nonlinear_shell_arclength.cpp. Clicking on a function +or class name will lead you to its reference documentation. + +\include nonlinear_shell_arclength.cpp + +*/ + +} \ No newline at end of file diff --git a/doc/nonlinear_shell_dynamic.dox b/doc/nonlinear_shell_dynamic.dox new file mode 100644 index 0000000..43b870c --- /dev/null +++ b/doc/nonlinear_shell_dynamic.dox @@ -0,0 +1,53 @@ +namespace gismo { +/** + +\page nonlinear_shell_dynamic Tutorial: Non-Linear dynamic analysis using Kirchhoff-Love shells + +This example shows how to perform a non-linear dynamic analysis using the \ref gsDynamicSolvers submodule of \ref gsStructuralAnalysis. The dynamic solvers provided by this submodule work on user-defined operators from \ref gsStructuralAnalysisOps, hence can be used with any element and in any application. In this example, we use the \ref gsKLShell. The goals of the tutorial are the following: +- Solve a non-linear quasi-static problem using \ref gsDynamicNewmark + +For detailed information about the construction of \ref gsThinShellAssembler and \ref gsMaterialMatrixBase, please consult the tutorials \ref linear_shell and \ref nonlinear_shell. + +\subsection nonlinear_shell_dynamic_Setup Problem set-up +The set-up of the problem to be solved for this example is identical to the one in \ref nonlinear_shell. The reader is referred to that example for a full overview of the code. + +\subsection nonlinear_shell_dynamic_Assembly Mass matrix assembly +The assembly of the linear stiffness matrix and the external force vector are similar as in \ref nonlinear_shell. To assemble the mass matrix, an additional assembly is performed using \c assembleMass(). The full assembly is given by: + +\snippet nonlinear_shell_dynamic.cpp Assemble linear part + +\subsection nonlinear_shell_dynamic_Jacobian Jacobian and Residual +To use the Newmark time integrator or any other time integrator from \ref gsDynamicSolvers, we need to define a Jacobian matrix and a residual, optionally time-independent. Both are defined from the \ref gsStructuralAnalysisOps, similar to \ref nonlinear_shell_static. + +In addition, an operator for mass and damping needs to be defined. These operators are simple lambda functions returning a constant matrix: + +\snippet nonlinear_shell_dynamic.cpp Define damping and mass matrices + +\subsection nonlinear_shell_dynamic_Solver Define the dynamic solver +Using the mass, damping and Jacobian matrices \f$M\f$, \f$C\f$ and \f$K_{NL}(\mathbf{u})\f$ and a residual vector \f$R(\mathbf{u})\f$, an non-linear dynamic solver can be defined. In the following, a \ref gsDynamicNewmark is defined: + +\snippet nonlinear_shell_dynamic.cpp Set dynamic solver + +The solution is initialized with zero-vectors, as follows: + +\snippet nonlinear_shell_dynamic.cpp Initialize solution + +\subsection nonlinear_shell_dynamic_Solve Solve the non-linear dynamic problem +The non-linear dynamic problem is solved by stepping over time. After some initializations, a simple time stepping loop is used: + +\snippet nonlinear_shell_dynamic.cpp Solve nonlinear problem + +It can be seen that the export of the solution field is similar as in \ref linear_shell \ref nonlinear_shell, but the Paraview files are added to a \ref gsParaviewCollection. This is saved at the end of the simulation as follows: + +\snippet nonlinear_shell_dynamic.cpp Save the paraview collection + +\section nonlinear_shell_dynamic_Annotated Annotated source file + +Here is the full file \c nonlinear_shell_dynamic.cpp. Clicking on a function +or class name will lead you to its reference documentation. + +\include nonlinear_shell_dynamic.cpp + +*/ + +} \ No newline at end of file diff --git a/doc/nonlinear_shell_static.dox b/doc/nonlinear_shell_static.dox new file mode 100644 index 0000000..7aa8857 --- /dev/null +++ b/doc/nonlinear_shell_static.dox @@ -0,0 +1,44 @@ +namespace gismo { +/** + +\page nonlinear_shell_static Tutorial: Non-Linear Kirchhoff-Love shell analysis using the gsStructuralAnalysis module + +This example solves a non-linear shell problem using the \ref gsThinShellAssembler and a static solver from \ref gsStructuralAnalysis. The goals of the tutorial are the following: +- Solve a geometrically non-linear problem using \ref gsStaticNewton + +For detailed information about the construction of \ref gsThinShellAssembler and \ref gsMaterialMatrixBase, please consult the tutorials \ref linear_shell and \ref nonlinear_shell. + +\subsection nonlinear_shell_static_Setup Problem set-up +The set-up of the problem to be solved for this example is identical to the one in \ref nonlinear_shell. The reader is referred to that example for a full overview of the code. + +\subsection nonlinear_shell_static_Jacobian Jacobian and Residual +In \ref nonlinear_shell, the Jacobian and Residual operators are already defined as \c std::function objects. In the present example, we use a similar approach, now defined through \ref gsStructuralAnalysisOps. In particular, the Jacobian and Residual operators are given by: + +\snippet nonlinear_shell_static.cpp Define nonlinear residual functions + +\subsection nonlinear_shell_static_Solver Define the static solver +The Newton-Raphson solver - \ref gsStaticNewton - used in this example is derived from \ref gsStaticBase. There are also other solvers available in the \ref gsStaticSolvers submodule. The Newton-Raphson solver needs to be defined using a linear stiffness matrix \f$K\f$, an external force vector \f$F\f$, a Jacobian matrix \f$K_{NL}(\mathbf{u})\f$ and a residual vector \f$R(\mathbf{u})\f$. The latter two are defined above in \ref nonlinear_shell_static_Jacobian, whereas the former two are defined like in \ref linear_shell_Assembler: + +\snippet nonlinear_shell_static.cpp Assemble linear part + +Consequently, the static solver is constructed as follows: + +\snippet nonlinear_shell_static.cpp Set static solver + +\subsection nonlinear_shell_static_Solve Solve the non-linear problem +The \ref gsStaticNewton simply performs Newton-Raphson iterations under the hood. Solving the non-linear system is simply done by using: + +\snippet nonlinear_shell_static.cpp Solve nonlinear problem + +Where the solution vector is obtained in the last line. As can be seen in the snipped above, the solver returns a status via \ref gsStatus, which can be used to check whether the solver converged or not. Evaluation and export of the solutions is similar as in \ref linear_shell and \ref nonlinear_shell. + +\section nonlinear_shell_static_Annotated Annotated source file + +Here is the full file \c nonlinear_shell_static.cpp. Clicking on a function +or class name will lead you to its reference documentation. + +\include nonlinear_shell_static.cpp + +*/ + +} \ No newline at end of file diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt new file mode 100644 index 0000000..689a545 --- /dev/null +++ b/examples/CMakeLists.txt @@ -0,0 +1,26 @@ +###################################################################### +## CMakeLists.txt --- gsStructuralAnalysis/examples +## This file is part of the G+Smo library. +## +## Author: Angelos Mantzaflaris, Hugo Verhelst +## Copyright (C) 2023 +###################################################################### + +# add example files +aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) +foreach(file ${FILES}) + add_gismo_executable(${file}) + get_filename_component(tarname ${file} NAME_WE) # name without extension + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + if(GISMO_BUILD_EXAMPLES) + set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + else(GISMO_BUILD_EXAMPLES) + set_target_properties(${tarname} PROPERTIES + FOLDER "${PROJECT_NAME}" + EXCLUDE_FROM_ALL TRUE) + endif(GISMO_BUILD_EXAMPLES) + add_dependencies(${PROJECT_NAME}-examples ${tarname}) + add_dependencies(${PROJECT_NAME}-all ${tarname}) + # install the example executables (optionally) + install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) +endforeach(file ${FILES}) \ No newline at end of file diff --git a/examples/example_DynamicBeamNL.cpp b/examples/example_DynamicBeamNL.cpp index fa2e8ca..fc546e6 100644 --- a/examples/example_DynamicBeamNL.cpp +++ b/examples/example_DynamicBeamNL.cpp @@ -18,7 +18,13 @@ #include #endif -#include +#include +#include +#include +#include +#include +#include + using namespace gismo; template @@ -96,21 +102,6 @@ int main (int argc, char** argv) gsInfo<<"EA:\t"<numDofs(); - gsMatrix<> uOld(N,1); - gsMatrix<> vOld(N,1); - gsMatrix<> aOld(N,1); - - uOld = assembler->projectL2(initialShape); + gsVector<> U(N); + gsVector<> V(N); + gsVector<> A(N); + // U.setZero(); + // V.setZero(); + // A.setZero(); + + U = assembler->projectL2(initialShape); // Set initial velocity to zero - vOld.setZero(); - aOld = assembler->projectL2(initialAcc); + V.setZero(); + A = assembler->projectL2(initialAcc); // Compute Error gsMultiPatch<> deformationIni = mp; - assembler->constructSolution(uOld,deformationIni); + assembler->constructSolution(U,deformationIni); deformationIni.patch(0).coefs() -= mp.patch(0).coefs(); gsField<> inifield(mp,deformationIni); gsInfo<<"Initial error:\t"<matrix(); // pre-assemble system assembler->assemble(); -// // set damping Matrix (same dimensions as M) -// C.setZero(M.rows(),M.cols()); +gsSparseMatrix<> C = gsSparseMatrix<>(assembler->numDofs(),assembler->numDofs()); +gsStructuralAnalysisOps::Damping_t Damping = [&C](const gsVector &, gsSparseMatrix & m) { m = C; return true; }; +gsStructuralAnalysisOps::Mass_t Mass = [&M]( gsSparseMatrix & m) { m = M; return true; }; + +gsDynamicBase * timeIntegrator; +if (method==1) + timeIntegrator = new gsDynamicExplicitEuler(Mass,Damping,Jacobian,Residual); +else if (method==2) + timeIntegrator = new gsDynamicImplicitEuler(Mass,Damping,Jacobian,Residual); +else if (method==3) + timeIntegrator = new gsDynamicNewmark(Mass,Damping,Jacobian,Residual); +else if (method==4) + timeIntegrator = new gsDynamicBathe(Mass,Damping,Jacobian,Residual); +else if (method==5) +{ + timeIntegrator = new gsDynamicWilson(Mass,Damping,Jacobian,Residual); + timeIntegrator->options().setReal("gamma",1.4); +} +else if (method==6) + timeIntegrator = new gsDynamicRK4(Mass,Damping,Jacobian,Residual); +else + GISMO_ERROR("Method "< timeIntegrator(M,Jacobian,Residual,dt); +timeIntegrator->options().setReal("DT",dt); +timeIntegrator->options().setReal("TolU",1e-3); +timeIntegrator->options().setSwitch("Verbose",true); -timeIntegrator.verbose(); -timeIntegrator.setTolerance(1e-6); -timeIntegrator.setMethod(methodName); -// if (quasiNewton) -// timeIntegrator.quasiNewton(); +timeIntegrator->setU(U); +timeIntegrator->setV(V); +timeIntegrator->setA(A); //------------------------------------------------------------------------------ // Initial Conditions //------------------------------------------------------------------------------ -gsMatrix<> uNew,vNew,aNew; -uNew = uOld; -vNew = vOld; -aNew = aOld; -timeIntegrator.setDisplacement(uNew); -timeIntegrator.setVelocity(vNew); -timeIntegrator.setAcceleration(aNew); - for (index_t i=0; istep(time,dt,U,V,A); + GISMO_ASSERT(status == gsStatus::Success,"Time integrator did not succeed"); - assembler->constructSolution(uNew,mp_def); + assembler->constructSolution(U,mp_def); + time = timeIntegrator->time(); - time = timeIntegrator.currentTime(); // Define manufactured solution char buffer_w_an[200]; sprintf(buffer_w_an,"%e*x/(24*%e)*(%e^3-2*%e*x^2+x^3)*cos(%e*pi*%e)",load,EI,length,length,omega,time); @@ -331,7 +332,7 @@ for (index_t i=0; i analytical("0","0",w_an,3); // Update solution and export - assembler->constructSolution(uNew,mp_def); + assembler->constructSolution(U,mp_def); mp_def.patch(0).coefs() -= mp.patch(0).coefs();// assuming 1 patch here gsField<> solField(mp,mp_def); std::string fileName = dirname + "/solution" + util::to_string(i); diff --git a/examples/example_DynamicShell.cpp b/examples/example_DynamicShell.cpp index 2ce9db8..d135a18 100644 --- a/examples/example_DynamicShell.cpp +++ b/examples/example_DynamicShell.cpp @@ -23,7 +23,7 @@ #include #endif -#include +#include #include using namespace gismo; @@ -101,10 +101,9 @@ int main (int argc, char** argv) gsInfo<<"Density:\t"< & result){re // C.setZero(M.rows(),M.cols()); // -gsTimeIntegrator timeIntegrator(M,K,TForce,dt); - -timeIntegrator.verbose(); -timeIntegrator.setTolerance(1e-6); -timeIntegrator.setMethod(methodName); +gsDynamicWilson timeIntegrator(Mass,Damping,Stiffness,TForce); +timeIntegrator.options().setReal("DT",dt); +timeIntegrator.options().setReal("TolU",1e-2); +timeIntegrator.options().setSwitch("Verbose",true); //------------------------------------------------------------------------------ // Initial Conditions @@ -269,9 +267,9 @@ gsMatrix<> uNew,vNew,aNew; uNew = uOld; vNew = vOld; aNew = aOld; -timeIntegrator.setDisplacement(uNew); -timeIntegrator.setVelocity(vNew); -timeIntegrator.setAcceleration(aNew); +timeIntegrator.setU(uNew); +timeIntegrator.setV(vNew); +timeIntegrator.setA(aNew); //------------------------------------------------------------------------------ // Nonlinear time integration @@ -280,12 +278,10 @@ real_t time; for (index_t i=0; i displacements = timeIntegrator.displacements(); + gsMatrix<> displacements = timeIntegrator.solutionU(); assembler->constructSolution(displacements,solution); @@ -302,7 +298,7 @@ for (index_t i=0; i res2; solution.patch(0).eval_into(v,res2); - time = timeIntegrator.currentTime(); + time = timeIntegrator.time(); std::ofstream file; file.open(wn,std::ofstream::out | std::ofstream::app); file << std::setprecision(10) diff --git a/examples/example_DynamicShellNL.cpp b/examples/example_DynamicShellNL.cpp index 033643e..bf72269 100644 --- a/examples/example_DynamicShellNL.cpp +++ b/examples/example_DynamicShellNL.cpp @@ -23,7 +23,11 @@ #include #endif -#include +#include +#include +#include +#include +#include using namespace gismo; // Choose among various shell examples, default = Thin Plate @@ -101,27 +105,12 @@ int main (int argc, char** argv) gsInfo<<"Density:\t"< Minv; - gsSparseMatrix<> M; + gsSparseMatrix<> M, C; gsSparseMatrix<> K; gsSparseMatrix<> K_T; @@ -269,19 +257,34 @@ gsStructuralAnalysisOps::TResidual_t Residual = [&assembler,&mp_def](gsM // Compute mass matrix (since it is constant over time) assembler->assembleMass(); M = assembler->matrix(); +C = gsSparseMatrix<>(assembler->numDofs(),assembler->numDofs()); // pre-assemble system assembler->assemble(); -// // set damping Matrix (same dimensions as M) -// C.setZero(M.rows(),M.cols()); - -gsTimeIntegrator timeIntegrator(M,Jacobian,Residual,dt); +gsStructuralAnalysisOps::Mass_t Mass = [&M]( gsSparseMatrix & m) { m = M; return true; }; +gsStructuralAnalysisOps::Damping_t Damping = [&C](const gsVector &, gsSparseMatrix & m) { m = C; return true; }; + +gsDynamicBase * timeIntegrator; +if (method==1) + timeIntegrator = new gsDynamicExplicitEuler(Mass,Damping,Jacobian,Residual); +else if (method==2) + timeIntegrator = new gsDynamicImplicitEuler(Mass,Damping,Jacobian,Residual); +else if (method==3) + timeIntegrator = new gsDynamicNewmark(Mass,Damping,Jacobian,Residual); +else if (method==4) + timeIntegrator = new gsDynamicBathe(Mass,Damping,Jacobian,Residual); +else if (method==5) +{ + timeIntegrator = new gsDynamicWilson(Mass,Damping,Jacobian,Residual); + timeIntegrator->options().setReal("gamma",1.4); +} +else + GISMO_ERROR("Method "<options().setReal("DT",dt); +timeIntegrator->options().setReal("TolU",1e-3); +timeIntegrator->options().setSwitch("Quasi",quasiNewton); +timeIntegrator->options().setSwitch("Verbose",true); //------------------------------------------------------------------------------ // Initial Conditions @@ -290,19 +293,18 @@ gsMatrix<> uNew,vNew,aNew; uNew = uOld; vNew = vOld; aNew = aOld; -timeIntegrator.setDisplacement(uNew); -timeIntegrator.setVelocity(vNew); -timeIntegrator.setAcceleration(aNew); +timeIntegrator->setU(uNew); +timeIntegrator->setV(vNew); +timeIntegrator->setA(aNew); real_t time; for (index_t i=0; istep(); if (status!=gsStatus::Success) GISMO_ERROR("Time integrator did not succeed"); - timeIntegrator.constructSolution(); - gsMatrix<> displacements = timeIntegrator.displacements(); + gsMatrix<> displacements = timeIntegrator->solutionU(); assembler->constructSolution(displacements,mp_def); @@ -319,7 +321,7 @@ for (index_t i=0; i res2; mp_def.patch(0).eval_into(v,res2); - time = timeIntegrator.currentTime(); + time = timeIntegrator->time(); std::ofstream file; file.open(wn,std::ofstream::out | std::ofstream::app); file << std::setprecision(10) @@ -337,6 +339,8 @@ collection.save(); delete materialMatrix; delete assembler; +delete timeIntegrator; + return result; } diff --git a/examples/example_DynamicShell_XBraid.cpp b/examples/example_DynamicShell_XBraid.cpp new file mode 100644 index 0000000..965aff0 --- /dev/null +++ b/examples/example_DynamicShell_XBraid.cpp @@ -0,0 +1,351 @@ +/** @file example_DynamicShell.cpp + + @brief Example for nonlinear time integration of a linear shell + + Fig 12 of: + Filho, L. A. D., & Awruch, A. M. (2004). + Geometrically nonlinear static and dynamic analysis of shells and plates using the eight-node hexahedral element with one-point quadrature. + Finite Elements in Analysis and Design. https://doi.org/10.1016/j.finel.2003.08.012 + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-..., TU Delft) +*/ + +#include + +#ifdef gsKLShell_ENABLED +#include +#include +#endif + +#include +#include +#include + +using namespace gismo; + +// template +// struct DynamicForce : public DynamicFunctor +// { +// DynamicForce(gsVector & F) +// : +// m_F(F) +// {} + +// gsVector m_F; + +// index_t rows() {return m_F.rows(); }; +// index_t cols() {return 1;}; + +// }; + + +// Choose among various shell examples, default = Thin Plate +#ifdef gsKLShell_ENABLED +int main (int argc, char** argv) +{ + // Input options + int numElevate = 1; + int numHref = 1; + bool plot = false; + + real_t thickness = 0.01576; + real_t width = 1; // Width of the strip is equal to 1. + real_t Area = thickness*width; + + real_t E_modulus = 1e7; + real_t PoissonRatio = 0.3; + real_t Density = 0.000245; + gsMultiPatch<> mp; + + real_t EA = E_modulus*Area; + real_t EI = 1.0/12.0*(width*math::pow(thickness,3))*E_modulus; + + int testCase = 0; + + int method = 2; // 1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe + + int result = 0; + std::string fn("surface/sphericalCap.xml"); + + bool write = false; + std::string wn; + + int steps = 100; + real_t tend = 3e-4; + + std::string assemberOptionsFile("options/solver_options.xml"); + + gsCmdLine cmd("Dynamics of a linear spherical cap."); + cmd.addString( "f", "file", "Input XML file for assembler options", assemberOptionsFile ); + cmd.addInt("r","hRefine", + "Number of dyadic h-refinement (bisection) steps to perform before solving", + numHref); + cmd.addInt("t", "testcase", + "Test case: 0: clamped-clamped, 1: pinned-pinned, 2: clamped-free", + testCase); + cmd.addInt("m", "method", + "1: Explicit Euler, 2: Implicit Euler, 3: Newmark, 4: Bathe", + method); + cmd.addInt("s", "steps", + "Number of time steps", + steps); + cmd.addReal("p", "endtime", + "End time of simulation", + tend); +// cmd.addString("g","geometry","File containing Geometry (.xml, .axl, .txt)",fn); + cmd.addInt("e","degreeElevation", + "Number of degree elevation steps to perform on the Geometry's basis before solving", + numElevate); + cmd.addSwitch("plot", "Plot result in ParaView format", plot); + cmd.addSwitch("write", "Write convergence data to file", write); + + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + gsFileData<> fd(assemberOptionsFile); + gsOptionList opts; + fd.getFirst(opts); + + gsInfo<<"Simulation time:\t"<(fn, mp); + + for(index_t i = 0; i< numElevate; ++i) + mp.patch(0).degreeElevate(); // Elevate the degree + + // h-refine + for(index_t i = 0; i< numHref; ++i) + mp.patch(0).uniformRefine(); + + gsInfo<<"Basis (patch 0): "<< mp.patch(0).basis() << "\n"; + +//------------------------------------------------------------------------------ +// Define Boundary conditions and Initial conditions +//------------------------------------------------------------------------------ + + gsBoundaryConditions<> BCs; + BCs.setGeoMap(mp); + gsPointLoads pLoads = gsPointLoads(); + + BCs.addCondition(boundary::north, condition_type::dirichlet, 0, 0, false, 0 ); // unknown 0 - x + BCs.addCondition(boundary::north, condition_type::clamped,0,0,false,2); + + BCs.addCondition(boundary::east, condition_type::dirichlet, 0, 0, false, 0 ); // unknown 0 - x + BCs.addCondition(boundary::east, condition_type::dirichlet, 0, 0, false, 1 ); // unknown 1 - y + BCs.addCondition(boundary::east, condition_type::dirichlet, 0, 0, false, 2 ); // unknown 2 - z + BCs.addCondition(boundary::east, condition_type::clamped,0,0,false,2); + + BCs.addCondition(boundary::south, condition_type::dirichlet, 0, 0, false, 1 ); // unknown 1 - y + BCs.addCondition(boundary::south, condition_type::clamped,0,0,false,2); + + BCs.addCondition(boundary::west, condition_type::clamped,0,0,false,2); + BCs.addCondition(boundary::west, condition_type::collapsed,0,0,false,2); + + gsVector<> point(2); point<< 0, 0 ; // Point where surface force is applied + gsVector<> loadvec (3); loadvec << 0, 0, -25 ; + pLoads.addLoad(point, loadvec, 0 ); + + gsVector<> tmp(3); + tmp<<0,0,0; + +//------------------------------------------------------------------------------ +// Define Beam Assembler and Assembler for L2-projection +//------------------------------------------------------------------------------ + + gsMultiPatch<> mp_def = mp, solution = mp; + gsMultiBasis<> dbasis(mp); + // Linear isotropic material model + gsConstantFunction<> force(tmp,3); + gsFunctionExpr<> t(std::to_string(thickness), 3); + gsFunctionExpr<> E(std::to_string(E_modulus),3); + gsFunctionExpr<> nu(std::to_string(PoissonRatio),3); + gsFunctionExpr<> rho(std::to_string(Density),3); + + std::vector*> parameters(2); + parameters[0] = &E; + parameters[1] = ν + + gsMaterialMatrixBase* materialMatrix; + + gsOptionList options; + options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",0); + options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",1); + materialMatrix = getMaterialMatrix<3,real_t>(mp,t,parameters,rho,options); + + gsThinShellAssemblerBase* assembler; + assembler = new gsThinShellAssembler<3, real_t, true >(mp,dbasis,BCs,force,materialMatrix); + + // Construct assembler object for dynamic computations + assembler->setOptions(opts); + assembler->setPointLoads(pLoads); + + assembler->assemble(); + size_t N = assembler->numDofs(); + gsMatrix<> uOld(N,1); + gsMatrix<> vOld(N,1); + gsMatrix<> aOld(N,1); + uOld.setZero(); + vOld.setZero(); + aOld.setZero(); + + +//------------------------------------------------------------------------------ +// Initiate mass and stiffness matrices and vectors for velocity, displacement and acceleration +//------------------------------------------------------------------------------ + + gsMatrix Minv; + gsSparseMatrix<> M; + gsSparseMatrix<> K; + gsSparseMatrix<> K_T; + gsVector<> F; + +//------------------------------------------------------------------------------ +// Nonlinear time integration +//------------------------------------------------------------------------------ +gsParaviewCollection collection(dirname + "/solution"); + +// Compute mass matrix (since it is constant over time) +assembler->assembleMass(); +M = assembler->matrix(); +// pre-assemble system +assembler->assemble(); +K = assembler->matrix(); +F = assembler->rhs(); +// Function for the Residual +gsStructuralAnalysisOps::Mass_t Mass; +gsStructuralAnalysisOps::Damping_t Damping; +gsStructuralAnalysisOps::Stiffness_t Stiffness; +gsStructuralAnalysisOps::TForce_t TForce; + +Mass = [&M]( gsSparseMatrix & result){result = M; return true;}; +Damping = [&M]( const gsVector & x, gsSparseMatrix & result){result = gsSparseMatrix(M.rows(),M.cols()); return true;}; +Stiffness = [&K]( gsSparseMatrix & result){result = K; return true;}; +TForce = [&F](real_t time, gsVector & result){result = F; return true;}; + +// // set damping Matrix (same dimensions as M) +// C.setZero(M.rows(),M.cols()); +// + +// DynamicForce myForce(...); + +gsDynamicNewmark timeIntegrator(Mass,Damping,Stiffness,TForce); // myForce(t) +timeIntegrator.options().setReal("DT",dt); +timeIntegrator.options().setReal("TolU",1e-2); +timeIntegrator.options().setSwitch("Verbose",false); + +gsMpiComm comm = gsMpi::init().worldComm(); +gsDynamicXBraid solver(&timeIntegrator,comm,0.0,tend,assembler->numDofs(),steps); + +//------------------------------------------------------------------------------ +// Initial Conditions +//------------------------------------------------------------------------------ +gsMatrix<> uNew,vNew,aNew; +uNew = uOld; +vNew = vOld; +aNew = aOld; +timeIntegrator.setU(uNew); +timeIntegrator.setV(vNew); +timeIntegrator.setA(aNew); + +//------------------------------------------------------------------------------ +// Nonlinear time integration +//------------------------------------------------------------------------------ + +solver.solve(); + + +// real_t time; +// for (index_t i=0; i displacements = timeIntegrator.solutionU(); + +// assembler->constructSolution(displacements,solution); + +// solution.patch(0).coefs() -= mp.patch(0).coefs();// assuming 1 patch here +// gsField<> solField(mp,solution); +// std::string fileName = dirname + "/solution" + util::to_string(i); +// gsWriteParaview<>(solField, fileName, 500); +// fileName = "solution" + util::to_string(i) + "0"; +// collection.addPart(fileName + ".vts",i); + +// if (write) +// { +// gsMatrix<> v(2,1); +// v<< 0.0,0.0; +// gsMatrix<> res2; +// solution.patch(0).eval_into(v,res2); +// time = timeIntegrator.time(); +// std::ofstream file; +// file.open(wn,std::ofstream::out | std::ofstream::app); +// file << std::setprecision(10) +// << time << "," +// << res2(2,0) <<"\n"; +// file.close(); +// } +// // Update solution with multipatch coefficients to generate geometry again +// solution.patch(0).coefs() += mp.patch(0).coefs();// assuming 1 patch here +// } +// collection.save(); + +delete materialMatrix; +delete assembler; + +return result; +} +#else//gsKLShell_ENABLED +int main(int argc, char *argv[]) +{ + gsWarn<<"G+Smo is not compiled with the gsKLShell module."; + return EXIT_FAILURE; +} +#endif \ No newline at end of file diff --git a/examples/example_PinchedMembrane_DWR.cpp b/examples/example_PinchedMembrane_DWR.cpp index f2fd8b9..56f6a19 100644 --- a/examples/example_PinchedMembrane_DWR.cpp +++ b/examples/example_PinchedMembrane_DWR.cpp @@ -286,9 +286,7 @@ int main(int argc, char *argv[]) gsThinShellAssemblerDWRBase * DWR; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - system(command); + gsFileManager::mkdir(dirname); gsParaviewCollection collection(dirname + "/" + "solution"); gsParaviewCollection errors(dirname + "/" + "error_elem_ref"); @@ -354,7 +352,6 @@ int main(int argc, char *argv[]) DWR->assemblePrimalL(); gsVector<> Force = DWR->primalL(); - typedef std::function (gsVector const &, real_t, gsVector const &) > ALResidual_t; // Function for the Jacobian gsStructuralAnalysisOps::Jacobian_t Jacobian = [&DWR,&mp_def](gsVector const &x, gsSparseMatrix & m) { diff --git a/examples/example_ShearWrinkling.cpp b/examples/example_ShearWrinkling.cpp index 3d3f5ed..e9aeab9 100644 --- a/examples/example_ShearWrinkling.cpp +++ b/examples/example_ShearWrinkling.cpp @@ -257,9 +257,7 @@ int main (int argc, char** argv) if (THB) dirname = dirname + "_THB"; - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - system(command); + gsFileManager::mkdir(dirname); // plot geometry if (plot) diff --git a/examples/gsElasticity_Modal.cpp b/examples/gsElasticity_Modal.cpp index 5fd1008..61a65b6 100644 --- a/examples/gsElasticity_Modal.cpp +++ b/examples/gsElasticity_Modal.cpp @@ -188,7 +188,7 @@ int main (int argc, char** argv) gsMassAssembler assemblerMass(mp,dbasis,BCs,g); assemblerMass.options() = assembler.options(); - assemblerMass.options().addReal("Density","Density of the material",1.); + assemblerMass.options().addReal("Density","Density of the material",Density); // Configure Structural Analsysis module assembler.assemble(); @@ -197,13 +197,20 @@ int main (int argc, char** argv) gsSparseMatrix<> M = assemblerMass.matrix(); gsModalSolver modal(K,M); - modal.verbose(); - + modal.options().setInt("solver",2); + modal.options().setInt("selectionRule",0); + modal.options().setInt("sortRule",4); + modal.options().setSwitch("verbose",true); + modal.options().setInt("ncvFac",2); + modal.options().setReal("shift",shift); + + gsStatus status; if (!sparse) - modal.compute(); + status = modal.compute(); else - modal.computeSparse(shift,10); + status = modal.computeSparse(10);//,2,Spectra::SortRule::LargestMagn,Spectra::SortRule::SmallestMagn); + GISMO_ENSURE(status == gsStatus::Success,"Buckling solver failed"); gsMatrix<> values = modal.values(); gsMatrix<> vectors = modal.vectors(); @@ -216,7 +223,7 @@ int main (int argc, char** argv) if (plot) { gsInfo<<"Plotting in Paraview...\n"; - system("mkdir -p ModalResults"); + gsFileManager::mkdir("ModalResults"); gsMatrix<> modeShape; gsMultiPatch<> displacement; gsParaviewCollection collection("ModalResults/modes_solid"); @@ -309,9 +316,9 @@ gsMultiPatch BrickDomain(int n, int m, int o, int p, int q ,int r, T L, T B, // Define a matrix with ones gsVector<> temp(len0); temp.setOnes(); - for (index_t l = 0; l < len2; l++) + for (size_t l = 0; l < len2; l++) { - for (index_t k = 0; k < len1; k++) + for (size_t k = 0; k < len1; k++) { index_t offset = l*len0*len1; // First column contains x-coordinates (length) diff --git a/examples/gsThinShell_ArcLength.cpp b/examples/gsThinShell_ArcLength.cpp index 3cc17d1..5463715 100644 --- a/examples/gsThinShell_ArcLength.cpp +++ b/examples/gsThinShell_ArcLength.cpp @@ -824,10 +824,7 @@ int main (int argc, char** argv) // SingularPoint = false; // } - std::string commands = "mkdir -p " + dirname; - const char *command = commands.c_str(); - int systemRet = system(command); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir(dirname); // plot geometry @@ -1146,9 +1143,9 @@ int main (int argc, char** argv) real_t S = Lold / 1e-3 / lambdas(0) / lambdas(2); real_t San = mu * (math::pow(lambdas(1),2)-1/lambdas(1)); - gsDebugVar(S); - gsDebugVar(San); - gsDebugVar(abs(S-San)); + gsInfo<<"S = "< deformation = solution; gsMatrix<> modeShape; @@ -588,8 +587,7 @@ int main (int argc, char** argv) if (write) { - int systemRet = system("mkdir -p ModalResults"); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); + gsFileManager::mkdir("ModalResults"); std::string wnM = "ModalResults/eigenvalues.txt"; writeToCSVfile(wnM,values); diff --git a/examples/snapping_element.cpp b/examples/snapping_element.cpp index b13838a..66ef776 100644 --- a/examples/snapping_element.cpp +++ b/examples/snapping_element.cpp @@ -360,7 +360,7 @@ int main(int argc, char *argv[]) if (arcLength->stabilityChange()) { gsInfo<<"Bifurcation spotted!"<<"\n"; - arcLength->computeSingularPoint(false); + arcLength->computeSingularPoint(Uold,Lold,false); arcLength->switchBranch(); dLb0 = dLb = dL; arcLength->setLength(dLb); diff --git a/examples/snapping_element3D.cpp b/examples/snapping_element3D.cpp index b96def7..11eb0a6 100644 --- a/examples/snapping_element3D.cpp +++ b/examples/snapping_element3D.cpp @@ -138,7 +138,7 @@ int main(int argc, char *argv[]) for (size_t p=0; p!=element.nPatches(); p++) { gsTensorBSpline<3,real_t>::uPtr tmp = gsNurbsCreator<>::lift3D(dynamic_cast &>(element.patch(p)),3e-3); - mp.addPatch(*tmp,element.patch(p).label()); + mp.addPatch(*tmp); } index_t topmid_ID = 14; mp.computeTopology(); @@ -381,7 +381,7 @@ else if (arcLength->stabilityChange()) { gsInfo<<"Bifurcation spotted!"<<"\n"; - arcLength->computeSingularPoint(false); + arcLength->computeSingularPoint(Uold,Lold,false); arcLength->switchBranch(); dLb0 = dLb = dL; arcLength->setLength(dLb); diff --git a/examples/snapping_example.cpp b/examples/snapping_example.cpp index cd25541..553776a 100644 --- a/examples/snapping_example.cpp +++ b/examples/snapping_example.cpp @@ -149,41 +149,50 @@ int main(int argc, char *argv[]) gsInfo<<"-----------------------------------Making the element------------------------------------------------\n"; gsInfo<<"-----------------------------------------------------------------------------------------------------\n"; gsMultiPatch<> element = makeElement(tw,tg,tb,ts,l,a,curves); - - gsWriteParaview(element,"element",1000,true); + auto elementlabels = element.getBoxProperty("label"); gsInfo<<"-----------------------------------------------------------------------------------------------------\n"; gsInfo<<"-----------------------------------Making the bottom block-------------------------------------------\n"; gsInfo<<"-----------------------------------------------------------------------------------------------------\n"; gsMultiPatch<> bottom = makeBottom(tw,tg,tb,ts,l,a,curves); - - gsWriteParaview(bottom,"bottom",1000,true); + auto bottomlabels = bottom.getBoxProperty("label"); gsInfo<<"-----------------------------------------------------------------------------------------------------\n"; gsInfo<<"-----------------------------------Making the top block-------------------------------------------\n"; gsInfo<<"-----------------------------------------------------------------------------------------------------\n"; gsMultiPatch<> top = makeTop(tw,tg,tb,ts,l,a,curves); - - gsWriteParaview(top,"top",1000,true); + auto toplabels = top.getBoxProperty("label"); gsInfo<<"-----------------------------------------------------------------------------------------------------\n"; gsInfo<<"-----------------------------------Making the geometry--------------------------------------------\n"; gsInfo<<"-----------------------------------------------------------------------------------------------------\n"; gsMultiPatch<> mp; + auto labels = mp.addBoxProperty("label",std::string()); + real_t dx = l; real_t dy = 2*tg+ts+tb; index_t topmid_ID; + index_t pIndex; + gsMultiPatch<> tmp; for (index_t kx = 0; kx!=Nx; kx++) { - gsMultiPatch<> tmp; tmp = bottom; gsNurbsCreator<>::shift2D(tmp,kx*dx,0); - mp.addPatches(tmp); + for (auto patch = tmp.begin(); patch != tmp.end(); patch++) + { + pIndex = mp.addPatch(**patch); + labels[pIndex] = bottomlabels[(*patch)->id()]; + } + for (index_t ky = 0; ky!=Ny; ky++) { tmp = element; gsNurbsCreator<>::shift2D(tmp,kx*dx,ky*dy); - mp.addPatches(tmp); + for (auto patch = tmp.begin(); patch != tmp.end(); patch++) + { + pIndex = mp.addPatch(**patch); + labels[pIndex] = elementlabels[(*patch)->id()]; + } } tmp = top; gsNurbsCreator<>::shift2D(tmp,kx*dx,Ny*dy); @@ -196,8 +205,11 @@ int main(int argc, char *argv[]) gsDebugVar(topmid_ID); } - mp.addPatches(tmp); - + for (auto patch = tmp.begin(); patch != tmp.end(); patch++) + { + pIndex = mp.addPatch(**patch); + labels[pIndex] = toplabels[(*patch)->id()]; + } } mp.computeTopology(); @@ -436,7 +448,7 @@ else if (arcLength->stabilityChange()) { gsInfo<<"Bifurcation spotted!"<<"\n"; - arcLength->computeSingularPoint(false); + arcLength->computeSingularPoint(Uold,Lold,false); arcLength->switchBranch(); dLb0 = dLb = dL; arcLength->setLength(dLb); diff --git a/examples/snapping_example3D.cpp b/examples/snapping_example3D.cpp index e045006..86a712e 100644 --- a/examples/snapping_example3D.cpp +++ b/examples/snapping_example3D.cpp @@ -184,7 +184,7 @@ int main(int argc, char *argv[]) for (size_t p = 0; p!=tmp.nPatches(); p ++) { gsTensorBSpline<3,real_t>::uPtr tmp2 = gsNurbsCreator<>::lift3D(dynamic_cast &>(tmp.patch(p)),b); - mp.addPatch(*tmp2,tmp.patch(p).label()); + mp.addPatch(*tmp2); } // mid @@ -195,7 +195,7 @@ int main(int argc, char *argv[]) for (size_t p = 0; p!=tmp.nPatches(); p ++) { gsTensorBSpline<3,real_t>::uPtr tmp2 = gsNurbsCreator<>::lift3D(dynamic_cast &>(tmp.patch(p)),b); - mp.addPatch(*tmp2,tmp.patch(p).label()); + mp.addPatch(*tmp2); } } @@ -212,7 +212,7 @@ int main(int argc, char *argv[]) for (size_t p = 0; p!=tmp.nPatches(); p ++) { gsTensorBSpline<3,real_t>::uPtr tmp2 = gsNurbsCreator<>::lift3D(dynamic_cast &>(tmp.patch(p)),b); - mp.addPatch(*tmp2,tmp.patch(p).label()); + mp.addPatch(*tmp2); } } mp.computeTopology(); @@ -458,7 +458,7 @@ else if (arcLength->stabilityChange()) { gsInfo<<"Bifurcation spotted!"<<"\n"; - arcLength->computeSingularPoint(false); + arcLength->computeSingularPoint(Uold,Lold,false); arcLength->switchBranch(); dLb0 = dLb = dL; arcLength->setLength(dLb); diff --git a/examples/snapping_example_shell.cpp b/examples/snapping_example_shell.cpp index 4a9c0de..f84aa73 100644 --- a/examples/snapping_example_shell.cpp +++ b/examples/snapping_example_shell.cpp @@ -473,16 +473,11 @@ int main(int argc, char *argv[]) gsInfo<<"Load step "<< k<<"\n"; gsStopwatch timer; gsStatus status = arcLength->step(); - if (status==gsStatus::Success) - gsDebug<<"Step successful\n"; - else if (status==gsStatus::NotConverged) - gsDebug<<"Not converged\n"; - else if (status==gsStatus::AssemblyError) - gsDebug<<"Assembly error\n"; - else if (status==gsStatus::SolverError) - gsDebug<<"Solver error\n"; - else if (status==gsStatus::OtherError) - gsDebug<<"Other error\n"; + if (status==gsStatus::Success) {gsDebug<<"Step successful\n";} + else if (status==gsStatus::NotConverged) {gsDebug<<"Not converged\n";} + else if (status==gsStatus::AssemblyError) {gsDebug<<"Assembly error\n";} + else if (status==gsStatus::SolverError) {gsDebug<<"Solver error\n";} + else if (status==gsStatus::OtherError) {gsDebug<<"Other error\n";} time += timer.stop(); if (status==gsStatus::NotConverged || status==gsStatus::AssemblyError) diff --git a/solvers/CMakeLists.txt b/solvers/CMakeLists.txt new file mode 100644 index 0000000..01fe917 --- /dev/null +++ b/solvers/CMakeLists.txt @@ -0,0 +1,26 @@ +###################################################################### +## CMakeLists.txt --- gsStructuralAnalysis/solvers +## This file is part of the G+Smo library. +## +## Author: Angelos Mantzaflaris, Hugo Verhelst +## Copyright (C) 2023 +###################################################################### + +# add solver files +aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) +foreach(file ${FILES}) + add_gismo_executable(${file}) + get_filename_component(tarname ${file} NAME_WE) # name without extension + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + if(GISMO_BUILD_EXAMPLES) + set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + else(GISMO_BUILD_EXAMPLES) + set_target_properties(${tarname} PROPERTIES + FOLDER "${PROJECT_NAME}" + EXCLUDE_FROM_ALL TRUE) + endif(GISMO_BUILD_EXAMPLES) + add_dependencies(${PROJECT_NAME}-solvers ${tarname}) + add_dependencies(${PROJECT_NAME}-all ${tarname}) + # install the example executables (optionally) + install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) +endforeach(file ${FILES}) \ No newline at end of file diff --git a/solvers/buckling_shell_XML.cpp b/solvers/buckling_shell_XML.cpp index 5723fd8..19e43a7 100644 --- a/solvers/buckling_shell_XML.cpp +++ b/solvers/buckling_shell_XML.cpp @@ -286,13 +286,12 @@ int main (int argc, char** argv) collection.save(); } - if (write) - { - int systemRet = system("mkdir -p BucklingResults"); - GISMO_ASSERT(systemRet!=-1,"Something went wrong with calling the system argument"); - std::string wnM = "BucklingResults/eigenvalues.txt"; - // writeToCSVfile(wnM,values); - } + // if (write) + // { + // gsFileManager::mkdir("BucklingResults"); + // std::string wnM = "BucklingResults/eigenvalues.txt"; + // writeToCSVfile(wnM,values); + // } return EXIT_SUCCESS; } #else//gsKLShell_ENABLED diff --git a/solvers/buckling_shell_multipatch_XML.cpp b/solvers/buckling_shell_multipatch_XML.cpp index 760d8e1..d4dd9f6 100644 --- a/solvers/buckling_shell_multipatch_XML.cpp +++ b/solvers/buckling_shell_multipatch_XML.cpp @@ -23,11 +23,13 @@ #include #endif +#ifdef gsUnstructuredSplines_ENABLED #include #include #include #include #include +#endif using namespace gismo; @@ -52,10 +54,10 @@ void writeToCSVfile(std::string name, gsMatrix<> matrix) } } -#ifdef gsKLShell_ENABLED -#ifdef gsUnstructuredSplines_ENABLED int main (int argc, char** argv) { +#ifdef gsKLShell_ENABLED +#ifdef gsUnstructuredSplines_ENABLED // Input options int numElevate = 0; int numRefine = 1; @@ -476,18 +478,12 @@ int main (int argc, char** argv) file.close(); } return EXIT_SUCCESS; -} #else//gsUnstructuredSplines_ENABLED -int main(int argc, char *argv[]) -{ gsWarn<<"G+Smo is not compiled with the gsUnstructuredSplines module."; return EXIT_FAILURE; -} #endif #else//gsKLShell_ENABLED -int main(int argc, char *argv[]) -{ gsWarn<<"G+Smo is not compiled with the gsKLShell module."; return EXIT_FAILURE; +#endif } -#endif \ No newline at end of file diff --git a/src/gsALMSolvers/gsALM_.cpp b/src/gsALMSolvers/gsALM_.cpp index 645a57a..e4c3905 100644 --- a/src/gsALMSolvers/gsALM_.cpp +++ b/src/gsALMSolvers/gsALM_.cpp @@ -17,6 +17,10 @@ namespace gismo { + + enum struct gsStatus; + template struct gsStructuralAnalysisOps; + CLASS_TEMPLATE_INST gsALMBase; CLASS_TEMPLATE_INST gsALMLoadControl; CLASS_TEMPLATE_INST gsALMRiks; diff --git a/src/gsALMSolvers/gsAPALM.hpp b/src/gsALMSolvers/gsAPALM.hpp index c27b4e2..f9dc657 100644 --- a/src/gsALMSolvers/gsAPALM.hpp +++ b/src/gsALMSolvers/gsAPALM.hpp @@ -298,7 +298,7 @@ gsAPALM::parallelSolve_impl() ID = std::get<0>(dataEntry); dataLevel = m_data.branch(branch).jobLevel(ID); bool success = m_data.branch(branch).getReferenceByID(ID,reference); - GISMO_ASSERT(success,"Reference not found"); + GISMO_ENSURE(success,"Reference not found"); this->_sendMainToWorker(m_workers.front(),stop); this->_sendMainToWorker(m_workers.front(), @@ -358,7 +358,7 @@ gsAPALM::parallelSolve_impl() ID = std::get<0>(dataEntry); dataLevel = m_data.branch(branch).jobLevel(ID); bool success = m_data.branch(branch).getReferenceByID(ID,reference); - GISMO_ASSERT(success,"Reference not found"); + GISMO_ENSURE(success,"Reference not found"); this->_sendMainToWorker(m_workers.front(),stop); this->_sendMainToWorker(m_workers.front(), branch, @@ -452,7 +452,7 @@ gsAPALM::parallelSolve_impl() ID = std::get<0>(dataEntry); dataLevel = m_data.branch(branch).jobLevel(ID); bool success = m_data.branch(branch).getReferenceByID(ID,reference); - GISMO_ASSERT(success,"Reference not found"); + GISMO_ENSURE(success,"Reference not found"); this->_correction(dataEntry, m_data.branch(branch).jobTimes(ID), dataLevel, @@ -600,7 +600,7 @@ gsAPALM::_solve_impl(index_t Nsteps) else { bool success = m_data.branch(branch).getReferenceByID(ID,reference); - GISMO_ASSERT(success,"Reference not found"); + GISMO_ENSURE(success,"Reference not found"); this->_sendMainToWorker(m_workers.front(),stop); this->_sendMainToWorker(m_workers.front(), @@ -720,7 +720,7 @@ gsAPALM::_solve_impl(index_t Nsteps) else { bool success = m_data.branch(branch).getReferenceByID(ID,reference); - GISMO_ASSERT(success,"Reference not found"); + GISMO_ENSURE(success,"Reference not found"); this->_sendMainToWorker(m_workers.front(),stop); this->_sendMainToWorker(m_workers.front(), @@ -910,7 +910,7 @@ gsAPALM::_solve_impl(index_t Nsteps) std::vector stepSolutions; T lowerDistance, upperDistance; bool success = m_data.branch(branch).getReferenceByID(ID,reference); - GISMO_ASSERT(success,"Reference not found"); + GISMO_ENSURE(success,"Reference not found"); this->_correction(dataEntry, m_data.branch(branch).jobTimes(ID), dataLevel, diff --git a/src/gsALMSolvers/gsAPALMData.hpp b/src/gsALMSolvers/gsAPALMData.hpp index 490869c..ba5976d 100644 --- a/src/gsALMSolvers/gsAPALMData.hpp +++ b/src/gsALMSolvers/gsAPALMData.hpp @@ -334,7 +334,7 @@ void gsAPALMData::submit(index_t ID, const std::vector & distan // check if the total distance matches the original time step T dt_tmp = std::accumulate(distances.begin(), std::prev(distances.end()), 0.0); - GISMO_ASSERT((Dt-dt_tmp)/Dt<1e-12,"Total distance of the computed intervals should be equal to the original interval length ("< +#include +#include + +namespace gismo +{ + +/** + @brief Performs the arc length method to solve a nonlinear system of equations. + + \tparam T coefficient type + + \ingroup gsStructuralAnalysis +*/ +template +class gsDynamicBase +{ +protected: + + typedef typename gsStructuralAnalysisOps::Force_t Force_t; + typedef typename gsStructuralAnalysisOps::TForce_t TForce_t; + typedef typename gsStructuralAnalysisOps::Residual_t Residual_t; + typedef typename gsStructuralAnalysisOps::TResidual_t TResidual_t; + typedef typename gsStructuralAnalysisOps::Mass_t Mass_t; + typedef typename gsStructuralAnalysisOps::TMass_t TMass_t; + typedef typename gsStructuralAnalysisOps::Damping_t Damping_t; + typedef typename gsStructuralAnalysisOps::TDamping_t TDamping_t; + typedef typename gsStructuralAnalysisOps::Stiffness_t Stiffness_t; + typedef typename gsStructuralAnalysisOps::Jacobian_t Jacobian_t; + typedef typename gsStructuralAnalysisOps::TJacobian_t TJacobian_t; + +public: + + virtual ~gsDynamicBase() {}; + + /// Constructor + gsDynamicBase( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const Force_t & Force + ) + : + m_mass(Mass), + m_damping(Damping), + m_stiffness(Stiffness), + m_force(Force) + { + m_Tmass = [this]( const T time, gsSparseMatrix & result) -> bool {return m_mass(result);}; + m_Tdamping = [this](gsVector const & x, const T time, gsSparseMatrix & result) -> bool {return m_damping(x,result);}; + m_Tjacobian = [this](gsVector const & x, const T time, gsSparseMatrix & result) -> bool {return m_stiffness(result);}; + m_Tforce = [this]( const T time, gsVector & result) -> bool {return m_force(result);}; + m_Tresidual = [ ](gsVector const & x, const T time, gsVector & result) -> bool {GISMO_ERROR("time-dependent residual not available");}; + _init(); + } + + /// Constructor + gsDynamicBase( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const TForce_t & TForce + ) + : + m_mass(Mass), + m_damping(Damping), + m_stiffness(Stiffness), + m_Tforce(TForce) + { + m_Tmass = [this]( const T time, gsSparseMatrix & result) -> bool {return m_mass(result);}; + m_Tdamping = [this](gsVector const & x, const T time, gsSparseMatrix & result) -> bool {return m_damping(x,result);}; + m_Tjacobian = [this](gsVector const & x, const T time, gsSparseMatrix & result) -> bool {return m_stiffness(result);}; + m_Tresidual = [ ](gsVector const & x, const T time, gsVector & result) -> bool {GISMO_ERROR("time-dependent residual not available");}; + _init(); + } + + /// Constructor + gsDynamicBase( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const Residual_t & Residual + ) + : + m_mass(Mass), + m_damping(Damping), + m_jacobian(Jacobian), + m_residual(Residual) + { + m_Tmass = [this]( const T time, gsSparseMatrix & result) -> bool {return m_mass(result);}; + m_Tdamping = [this](gsVector const & x, const T time, gsSparseMatrix & result) -> bool {return m_damping(x,result);}; + m_Tjacobian = [this](gsVector const & x, const T time, gsSparseMatrix & result) -> bool {return m_jacobian(x,result);}; + m_Tforce = [ ]( const T time, gsVector & result) -> bool {GISMO_ERROR("time-dependent force not available");}; + m_Tresidual = [this](gsVector const & x, const T time, gsVector & result) -> bool {return m_residual(x,result);}; + _init(); + } + + /// Constructor + gsDynamicBase( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const TResidual_t & TResidual + ) + : + m_mass(Mass), + m_damping(Damping), + m_jacobian(Jacobian), + m_Tresidual(TResidual) + { + m_Tmass = [this]( const T time, gsSparseMatrix & result) -> bool {return m_mass(result);}; + m_Tjacobian = [this](gsVector const & x, const T time, gsSparseMatrix & result) -> bool {return m_jacobian(x,result);}; + m_Tdamping = [this](gsVector const & x, const T time, gsSparseMatrix & result) -> bool {return m_damping(x,result);}; + _init(); + } + + /// Constructor + gsDynamicBase( + const Mass_t & Mass, + const Damping_t & Damping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + m_mass(Mass), + m_damping(Damping), + m_Tjacobian(TJacobian), + m_Tresidual(TResidual) + { + m_Tmass = [this]( const T time, gsSparseMatrix & result) -> bool {return m_mass(result);}; + m_Tdamping = [this](gsVector const & x, const T time, gsSparseMatrix & result) -> bool {return m_damping(x,result);}; + _init(); + } + + /// Constructor + gsDynamicBase( + const TMass_t & TMass, + const TDamping_t & TDamping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + m_Tmass(TMass), + m_Tdamping(TDamping), + m_Tjacobian(TJacobian), + m_Tresidual(TResidual) + { + _init(); + } + + gsDynamicBase() {_init();} + +protected: + void _init() + { + m_time = 0; + // initialize variables + m_numIterations = -1; + defaultOptions(); + + m_status = gsStatus::NotStarted; + } + +// General functions +public: + + virtual gsStatus status() { return m_status; } + + // Returns the number of DoFs + // virtual index_t numDofs() {return m_forcing.size();} + + // Returns the current time step + virtual T getTimeStep() const {return m_options.getReal("DT"); } + + /// Perform one arc-length step + virtual gsStatus step(T dt) + { + gsStatus status = this->_step(m_time,dt,m_U,m_V,m_A); + m_time += dt; + return status; + } + + virtual gsStatus step() + { + return this->step(m_options.getReal("DT")); + } + + virtual gsStatus step(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const + { + return _step(t,dt,U,V,A); + } + + /// Set time step to \a dt + virtual void setTimeStep(T dt) + { + m_options.setReal("DT",dt); + } + + // Output + /// True if the Arc Length method converged + virtual bool converged() const {return m_status==gsStatus::Success;} + + /// Returns the number of Newton iterations performed + virtual index_t numIterations() const { return m_numIterations;} + + virtual const gsVector & displacements() const {return this->solutionU();} + virtual const gsVector & velocities() const {return this->solutionV();} + virtual const gsVector & accelerations() const {return this->solutionA();} + + virtual const gsVector & solutionU() const {return m_U;} + virtual const gsVector & solutionV() const {return m_V;} + virtual const gsVector & solutionA() const {return m_A;} + + virtual const T & time() {return m_time;} + + virtual void setU(const gsVector & U) {m_U = U;} + virtual void setV(const gsVector & V) {m_V = V;} + virtual void setA(const gsVector & A) {m_A = A;} + + virtual void setDisplacements(const gsVector & U) {this->setU(U);} + virtual void setVelocities(const gsVector & V) {this->setV(V);} + virtual void setAccelerations(const gsVector & A) {this->setA(A);} + + /// Access the options + virtual gsOptionList & options() {return m_options;}; + + /// Set the options to \a options + virtual void setOptions(gsOptionList options) {m_options.update(options,gsOptionList::addIfUnknown); }; + + /// Return the options into \a options + virtual const void options_into(gsOptionList options) {options = m_options;}; + + /// Return the number of degrees of freedom + virtual index_t numDofs() {return m_numDofs; } +// ------------------------------------------------------------------------------------------------------------ +// ---------------------------------------Computations--------------------------------------------------------- +// ------------------------------------------------------------------------------------------------------------ +protected: + /// Set default options + virtual void defaultOptions(); + + /// Compute the residual + virtual void _computeForce(const T time, gsVector & F) const + { + if (!m_Tforce(time,F)) + throw 2; + } + + /// Compute the residual + virtual void _computeResidual(const gsVector & U, const T time, gsVector & R) const + { + if (!m_Tresidual(U,time,R)) + throw 2; + } + + /// Compute the mass matrix + virtual void _computeMass(const T time, gsSparseMatrix & M) const + { + if (!m_Tmass(time,M)) + throw 2; + } + + /// Compute the mass matrix + virtual void _computeMassInverse(const gsSparseMatrix & M, gsSparseMatrix & Minv) const + { + if ((m_mass==nullptr) || (m_massInv.rows()==0 || m_massInv.cols()==0)) // then mass is time-dependent or the mass inverse is not stored, compute it + { + gsSparseMatrix eye(M.rows(), M.cols()); + eye.setIdentity(); + gsSparseSolver<>::LU solver(M); + gsMatrix MinvI = solver.solve(eye); + m_massInv = Minv = MinvI.sparseView(); + } + else + Minv = m_massInv; + } + + /// Compute the damping matrix + virtual void _computeDamping(const gsVector & U, const T time, gsSparseMatrix & C) const + { + if (!m_Tdamping(U,time,C)) + throw 2; + } + + /// Compute the Jacobian matrix + virtual void _computeJacobian(const gsVector & U, const T time, gsSparseMatrix & K) const + { + if (!m_Tjacobian(U,time,K)) + throw 2; + } + +// Purely virtual functions +protected: + /// Initialize the ALM + virtual gsStatus _step(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const = 0; + +protected: + + // Number of degrees of freedom + index_t m_numDofs; + + Mass_t m_mass; + mutable gsSparseMatrix m_massInv; + TMass_t m_Tmass; + + Damping_t m_damping; + TDamping_t m_Tdamping; + + Stiffness_t m_stiffness; + + Jacobian_t m_jacobian; + TJacobian_t m_Tjacobian; + + Force_t m_force; + TForce_t m_Tforce; + + Residual_t m_residual; + TResidual_t m_Tresidual; + + mutable typename gsSparseSolver::uPtr m_solver; // Cholesky by default + +protected: + + mutable gsOptionList m_options; + +protected: + + + /// Number of iterations performed + mutable index_t m_numIterations; + + // /// Maximum number of Arc Length iterations allowed + // index_t m_maxIterations; + + // /// Number of desired iterations + // index_t m_desiredIterations; + + // /// Time step + // T m_dt; + + /// Time + T m_time; + + // /// Tolerance value to decide convergence + // T m_tolerance; + + // /// Tolerance value to decide convergence - Force criterion + // T m_toleranceF; + + // /// Tolerance value to decide convergence - Displacement criterion + // T m_toleranceU; + + // bool m_verbose; + bool m_initialized; + + // bool m_quasiNewton; + // index_t m_quasiNewtonInterval; + + gsStatus m_status; + +protected: + + // /// Current update + // gsVector m_DeltaU; + // gsVector m_DeltaV; + // gsVector m_DeltaA; + + // /// Iterative update + // gsVector m_deltaU; + // gsVector m_deltaV; + // gsVector m_deltaA; + + // Previous update + gsVector m_DeltaUold; + gsVector m_DeltaVold; + gsVector m_DeltaAold; + + /// Displacement vector (present, at previously converged point) + gsVector m_U, m_Uold; + gsVector m_V, m_Vold; + gsVector m_A, m_Aold; + + /// Time + T m_t, m_tprev; +}; + + +} // namespace gismo + +#ifndef GISMO_BUILD_LIB +#include GISMO_HPP_HEADER(gsDynamicBase.hpp) +#endif diff --git a/src/gsDynamicSolvers/gsDynamicBase.hpp b/src/gsDynamicSolvers/gsDynamicBase.hpp new file mode 100644 index 0000000..8924231 --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicBase.hpp @@ -0,0 +1,36 @@ +/** @file gsDynamicBase.hpp + + @brief Base class to perform time integration of second-order structural dynamics systems + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-..., TU Delft) +*/ + +#pragma once + +namespace gismo +{ + +template +void gsDynamicBase::defaultOptions() +{ + m_options.addInt ("MaxIter","Maximum iterations",100); + m_options.addReal("TolF","Tolerance",1e-3); + m_options.addReal("TolU","Tolerance",1e-6); + + m_options.addReal("DT","Time step",1e-2); + + m_options.addSwitch("Quasi","Use Quasi Newton method",false); + m_options.addInt ("QuasiIterations","Number of iterations for quasi newton method",1); + + m_options.addString("Solver","Sparse linear solver", "SimplicialLDLT"); + + m_options.addSwitch ("Verbose","Verbose output",false); +} + +} // namespace gismo diff --git a/src/gsDynamicSolvers/gsDynamicBathe.h b/src/gsDynamicSolvers/gsDynamicBathe.h new file mode 100644 index 0000000..3e63d7f --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicBathe.h @@ -0,0 +1,184 @@ + /** @file gsDynamicBathe.h + + @brief Class to perform time integration of second-order structural dynamics systems using the Explicit Euler method + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-..., TU Delft) + + TODO (June 2023): + * Change inputs to const references! +*/ + +#pragma once +#include + +#include +#include +#include + +namespace gismo +{ + +/** + @brief Performs the arc length method to solve a nonlinear system of equations. + + \tparam T coefficient type + + \ingroup gsStructuralAnalysis +*/ +template +class gsDynamicBathe : public gsDynamicNewmark +{ + typedef gsDynamicNewmark Base; + +protected: + + typedef typename gsStructuralAnalysisOps::Force_t Force_t; + typedef typename gsStructuralAnalysisOps::TForce_t TForce_t; + typedef typename gsStructuralAnalysisOps::Residual_t Residual_t; + typedef typename gsStructuralAnalysisOps::TResidual_t TResidual_t; + typedef typename gsStructuralAnalysisOps::Mass_t Mass_t; + typedef typename gsStructuralAnalysisOps::TMass_t TMass_t; + typedef typename gsStructuralAnalysisOps::Damping_t Damping_t; + typedef typename gsStructuralAnalysisOps::TDamping_t TDamping_t; + typedef typename gsStructuralAnalysisOps::Stiffness_t Stiffness_t; + typedef typename gsStructuralAnalysisOps::Jacobian_t Jacobian_t; + typedef typename gsStructuralAnalysisOps::TJacobian_t TJacobian_t; + +public: + + virtual ~gsDynamicBathe() {}; + + /// Constructor + gsDynamicBathe( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const Force_t & Force + ) + : + Base(Mass,Damping,Stiffness,Force) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicBathe( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const TForce_t & TForce + ) + : + Base(Mass,Damping,Stiffness,TForce) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicBathe( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const Residual_t & Residual + ) + : + Base(Mass,Damping,Jacobian,Residual) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicBathe( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,Jacobian,TResidual) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicBathe( + const Mass_t & Mass, + const Damping_t & Damping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,TJacobian,TResidual) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicBathe( + const TMass_t & TMass, + const TDamping_t & TDamping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(TMass,TDamping,TJacobian,TResidual) + { + this->defaultOptions(); + } + + /// Set default options + virtual void defaultOptions() override; + +// General functions +protected: + + gsStatus _step(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const override + { + gsStatus status = gsStatus::NotStarted; + status = _step_impl<_NL>(t,dt,U,V,A); + return status; + } + + void _initOutput() const; + void _stageOutput(index_t stage) const; + void _stepOutput(const index_t it, const T resnorm, const T updatenorm) const; + + gsSparseMatrix m_sysmat; + + using Base::_computeForce; + using Base::_computeResidual; + using Base::_computeMass; + using Base::_computeMassInverse; + using Base::_computeDamping; + using Base::_computeJacobian; + +protected: + + using Base::m_solver; + + using Base::m_numIterations; + + using Base::m_options; + + +private: + template + typename std::enable_if<(_nonlinear==false), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; + + template + typename std::enable_if<(_nonlinear==true), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; +}; + +} // namespace gismo + +#ifndef GISMO_BUILD_LIB +#include GISMO_HPP_HEADER(gsDynamicBathe.hpp) +#endif diff --git a/src/gsDynamicSolvers/gsDynamicBathe.hpp b/src/gsDynamicSolvers/gsDynamicBathe.hpp new file mode 100644 index 0000000..6675818 --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicBathe.hpp @@ -0,0 +1,211 @@ +/** @file gsDynamicBathe.hpp + + @brief Performs the arc length method to solve a nonlinear equation system. + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-..., TU Delft) +*/ + +#pragma once + +namespace gismo +{ + +template +void gsDynamicBathe::defaultOptions() +{ + Base::defaultOptions(); + m_options.addReal("alpha","Beta parameter for Bathe's method, such that 0 =< 2 beta =< 1",0.25); + m_options.addReal("delta","Beta parameter for Bathe's method, such that 0 =< delta =< 1",0.50); + m_options.addReal("gamma","Beta parameter for Bathe's method, such that 0 =< gamma =< 1",0.50); +} + + +template +template +typename std::enable_if<(_nonlinear==false), gsStatus>::type +gsDynamicBathe::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + T gamma = m_options.getReal("gamma"); + + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + gsVector Ustep = U; + gsVector Vstep = V; + gsVector Astep = A; + + /// Stage 1: A Newmark step with DT=gamma*dt + this->_stageOutput(1); + Base::_step(t,gamma*dt,Ustep,Vstep,Astep); + + /// Stage 2: A Newmark step with DT=gamma*dt + gsVector F; + gsSparseMatrix M, C, K; + + // Computed at t=t0+dt + this->_computeMass(t+dt,M); + this->_computeForce(t+dt,F); + this->_computeDamping(U,t+dt,C); + this->_computeJacobian(U,t+dt,K); + + T c1 = (1-gamma)/(gamma*dt); + T c2 = (-1)/((1-gamma)*gamma*dt); + T c3 = (2-gamma)/((1-gamma)*dt); + + gsSparseMatrix lhs = K + c3*c3*M + c3*C; + gsMatrix rhs = F - M*(c1*Vold+c2*Vstep+c1*c3*Uold+c3*c2*Ustep) - C*(c2*Ustep+c1*Uold); + + this->_stageOutput(2); + this->_initOutput(); + typename gsSparseSolver::uPtr solver; + solver = gsSparseSolver::get( m_options.getString("Solver") ); + solver->compute(lhs); + U = solver->solve(rhs); + V = c1*Uold + c2*Ustep + c3*U; + A = c1*Vold + c2*Vstep + c3*V; + + this->_stepOutput(0,(F - K*U - M * A - C * V).norm(),U.norm()); + if (math::isinf(U.norm()) || math::isnan(U.norm())) + return gsStatus::NotConverged; + else + return gsStatus::Success; +} + + +template +template +typename std::enable_if<(_nonlinear==true), gsStatus>::type +gsDynamicBathe::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + T gamma = m_options.getReal("gamma"); + + T c1 = (1-gamma)/(gamma*dt); + T c2 = (-1)/((1-gamma)*gamma*dt); + T c3 = (2-gamma)/((1-gamma)*dt); + + + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + gsVector Ustep = U; + gsVector Vstep = V; + gsVector Astep = A; + + /// Stage 1: A Newmark step with DT=gamma*dt + this->_stageOutput(1); + Base::_step(t,gamma*dt,Ustep,Vstep,Astep); + + /// Stage 2: A Newmark step with DT=gamma*dt + gsVector R; + gsSparseMatrix M, C, K; + gsSparseMatrix lhs; + gsMatrix rhs; + gsMatrix dU; + + // Computed at t=t0+dt + this->_computeMass(t+dt,M); + this->_computeDamping(U,t+dt,C); + this->_computeResidual(U,t+dt,R); + + U = Ustep; + V = Vstep; + A = Astep; + + rhs = R - M*(c1*Vold+c2*Vstep+c1*c3*Uold+c3*c2*Ustep+c3*c3*U) - C*(c1*Uold+c2*Ustep+c3*U); + + T tolU = m_options.getReal("TolU"); + T tolF = m_options.getReal("TolF"); + T updateNorm = 10.0*tolU; + T residualNorm = rhs.norm(); + T residualNorm0 = (residualNorm!=0) ? residualNorm : 1; + this->_initOutput(); + T Unorm, dUnorm; + for (index_t numIterations = 0; numIterations < m_options.getInt("MaxIter"); ++numIterations) + { + if ((!m_options.getSwitch("Quasi")) || ((numIterations==0) || (numIterations % m_options.getInt("QuasiIterations") == 0))) + { + // Computed at t=t0+dt + this->_computeDamping(U,t+dt,C); + this->_computeJacobian(U,t+dt,K); + } + + lhs = K + c3*c3*M + c3*C; + + typename gsSparseSolver::uPtr solver; + solver = gsSparseSolver::get( m_options.getString("Solver") ); + solver->compute(lhs); + // U = solver->solve(rhs); + // V = delta/( dt*alpha ) * (U-Uold) - Vold; + // A = 1/( dt*dt*alpha ) * (U-Uold-Vold*dt) - Aold; + + dU = solver->solve(rhs); + U += dU; + + Unorm = U.norm(); + dUnorm = dU.norm(); + updateNorm = (Unorm!=0) ? dUnorm/Unorm : dUnorm; + + this->_computeResidual(U,t+dt,R); + rhs = R - M*(c1*Vold+c2*Vstep+c1*c3*Uold+c3*c2*Ustep+c3*c3*U) - C*(c1*Uold+c2*Ustep+c3*U); + residualNorm = rhs.norm() / residualNorm0; + updateNorm = dU.norm() / U.norm(); + + this->_stepOutput(numIterations,residualNorm,updateNorm); + + if ( (updateNorm +void gsDynamicBathe::_stageOutput(index_t stage) const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"Stage "< +void gsDynamicBathe::_initOutput() const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< +void gsDynamicBathe::_stepOutput(const index_t it, const T resnorm, const T updatenorm) const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< +#include +#include + +namespace gismo +{ + +/** + @brief Performs the arc length method to solve a nonlinear system of equations. + + \tparam T coefficient type + + \ingroup gsStructuralAnalysis +*/ +template +class gsDynamicExplicitEuler : public gsDynamicBase +{ + typedef gsDynamicBase Base; + +protected: + + typedef typename gsStructuralAnalysisOps::Force_t Force_t; + typedef typename gsStructuralAnalysisOps::TForce_t TForce_t; + typedef typename gsStructuralAnalysisOps::Residual_t Residual_t; + typedef typename gsStructuralAnalysisOps::TResidual_t TResidual_t; + typedef typename gsStructuralAnalysisOps::Mass_t Mass_t; + typedef typename gsStructuralAnalysisOps::TMass_t TMass_t; + typedef typename gsStructuralAnalysisOps::Damping_t Damping_t; + typedef typename gsStructuralAnalysisOps::TDamping_t TDamping_t; + typedef typename gsStructuralAnalysisOps::Stiffness_t Stiffness_t; + typedef typename gsStructuralAnalysisOps::Jacobian_t Jacobian_t; + typedef typename gsStructuralAnalysisOps::TJacobian_t TJacobian_t; + +public: + + virtual ~gsDynamicExplicitEuler() {}; + + /// Constructor + gsDynamicExplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const Force_t & Force + ) + : + Base(Mass,Damping,Stiffness,Force) + {} + + /// Constructor + gsDynamicExplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const TForce_t & TForce + ) + : + Base(Mass,Damping,Stiffness,TForce) + {} + + /// Constructor + gsDynamicExplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const Residual_t & Residual + ) + : + Base(Mass,Damping,Jacobian,Residual) + {} + + /// Constructor + gsDynamicExplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,Jacobian,TResidual) + {} + + /// Constructor + gsDynamicExplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,TJacobian,TResidual) + {} + + /// Constructor + gsDynamicExplicitEuler( + const TMass_t & TMass, + const TDamping_t & TDamping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(TMass,TDamping,TJacobian,TResidual) + {} + +// General functions +protected: + + gsStatus _step(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const override + { + gsStatus status = gsStatus::NotStarted; + status = _step_impl<_NL>(t,dt,U,V,A); + return status; + } + + void _initOutput() const; + void _stepOutput(const index_t it, const T resnorm, const T updatenorm) const; + + gsSparseMatrix m_sysmat; + + using Base::_computeForce; + using Base::_computeResidual; + using Base::_computeMass; + using Base::_computeMassInverse; + using Base::_computeDamping; + using Base::_computeJacobian; + +protected: + + using Base::m_solver; + + using Base::m_numIterations; + + using Base::m_options; + + +private: + template + typename std::enable_if<(_nonlinear==false), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; + + template + typename std::enable_if<(_nonlinear==true), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; +}; + +} // namespace gismo + +#ifndef GISMO_BUILD_LIB +#include GISMO_HPP_HEADER(gsDynamicExplicitEuler.hpp) +#endif diff --git a/src/gsDynamicSolvers/gsDynamicExplicitEuler.hpp b/src/gsDynamicSolvers/gsDynamicExplicitEuler.hpp new file mode 100644 index 0000000..8ddedc3 --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicExplicitEuler.hpp @@ -0,0 +1,151 @@ +/** @file gsDynamicExplicitEuler.hpp + + @brief Performs the arc length method to solve a nonlinear equation system. + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-..., TU Delft) +*/ + +#pragma once + +namespace gismo +{ + +// template +// gsStatus gsDynamicExplicitEuler::_step(const T dt) +// { +// gsVector R; +// gsSparseMatrix M, C; +// _computeResidual(m_U,m_time,R); +// _computeDamping(m_U,m_time,C); +// _computeMass(m_time,M); + +// m_Uold = m_U; +// m_Vold = m_V; +// m_Aold = m_A; + +// m_U = m_Uold + m_dt * m_Vold; +// // if (m_lumped) +// // { +// // gsVector M = _computeMassLumped(m_t); +// // m_DeltaV = m_dt * M.cwiseProduct(R - C * m_Vold);// element-wise; +// // } +// // else +// // { +// // gsSparseMatrix M = _computeMass(m_t); +// m_solver->compute(M); +// m_V = m_Vold + m_dt * m_solver->solve( R - C * m_Vold ); +// // } + +// m_time += m_dt; + +// return gsStatus::Success; +// } + +template +template +typename std::enable_if<(_nonlinear==false), gsStatus>::type +gsDynamicExplicitEuler::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + index_t N = U.rows(); + gsVector sol(2*N); + sol.topRows(N) = Uold; + sol.bottomRows(N) = Vold; + + gsVector F; + gsSparseMatrix M, Minv, C, K; + + // Computed at t=t0 + this->_computeMass(t,M); + this->_computeMassInverse(M,Minv); + this->_computeForce(t,F); + this->_computeDamping(U,t,C); + this->_computeJacobian(U,t,K); + + this->_initOutput(); + sol.topRows(N) += dt * Vold; + sol.bottomRows(N) += dt * Minv * (F - K * Uold - C * Vold); + this->_stepOutput(0,sol.norm(),0.); + gsDebugVar(sol.transpose()); + + U = sol.topRows(N); + V = sol.bottomRows(N); + + if (math::isinf(sol.norm()) || math::isnan(sol.norm())) + return gsStatus::NotConverged; + else + return gsStatus::Success; +} + + +template +template +typename std::enable_if<(_nonlinear==true), gsStatus>::type +gsDynamicExplicitEuler::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + index_t N = U.rows(); + gsVector sol(2*N); + sol.topRows(N) = Uold; + sol.bottomRows(N) = Vold; + + gsVector R; + gsSparseMatrix M, Minv, C, K; + + // Computed at t=t0 + this->_computeMass(t,M); + this->_computeMassInverse(M,Minv); + this->_computeDamping(Uold,t,C); + this->_computeResidual(Uold,t,R); + + this->_initOutput(); + sol.topRows(N) += dt * Vold; + sol.bottomRows(N) += dt * Minv * ( - R - C * Vold); + this->_stepOutput(0,sol.norm(),0.); + + U = sol.topRows(N); + V = sol.bottomRows(N); + + if (math::isinf(sol.norm()) || math::isnan(sol.norm())) + return gsStatus::NotConverged; + else + return gsStatus::Success; +} + +template +void gsDynamicExplicitEuler::_initOutput() const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< +void gsDynamicExplicitEuler::_stepOutput(const index_t it, const T resnorm, const T updatenorm) const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< +#include +#include + +namespace gismo +{ + +/** + @brief Performs the arc length method to solve a nonlinear system of equations. + + \tparam T coefficient type + + \ingroup gsStructuralAnalysis +*/ +template +class gsDynamicImplicitEuler : public gsDynamicBase +{ + typedef gsDynamicBase Base; + +protected: + + typedef typename gsStructuralAnalysisOps::Force_t Force_t; + typedef typename gsStructuralAnalysisOps::TForce_t TForce_t; + typedef typename gsStructuralAnalysisOps::Residual_t Residual_t; + typedef typename gsStructuralAnalysisOps::TResidual_t TResidual_t; + typedef typename gsStructuralAnalysisOps::Mass_t Mass_t; + typedef typename gsStructuralAnalysisOps::TMass_t TMass_t; + typedef typename gsStructuralAnalysisOps::Damping_t Damping_t; + typedef typename gsStructuralAnalysisOps::TDamping_t TDamping_t; + typedef typename gsStructuralAnalysisOps::Stiffness_t Stiffness_t; + typedef typename gsStructuralAnalysisOps::Jacobian_t Jacobian_t; + typedef typename gsStructuralAnalysisOps::TJacobian_t TJacobian_t; + +public: + + virtual ~gsDynamicImplicitEuler() {}; + + /// Constructor + gsDynamicImplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const Force_t & Force + ) + : + Base(Mass,Damping,Stiffness,Force) + {} + + /// Constructor + gsDynamicImplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const TForce_t & TForce + ) + : + Base(Mass,Damping,Stiffness,TForce) + {} + + /// Constructor + gsDynamicImplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const Residual_t & Residual + ) + : + Base(Mass,Damping,Jacobian,Residual) + {} + + /// Constructor + gsDynamicImplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,Jacobian,TResidual) + {} + + /// Constructor + gsDynamicImplicitEuler( + const Mass_t & Mass, + const Damping_t & Damping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,TJacobian,TResidual) + {} + + /// Constructor + gsDynamicImplicitEuler( + const TMass_t & TMass, + const TDamping_t & TDamping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(TMass,TDamping,TJacobian,TResidual) + {} + +// General functions +protected: + + gsStatus _step(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const override + { + gsStatus status = gsStatus::NotStarted; + status = _step_impl<_NL>(t,dt,U,V,A); + return status; + } + + void _initOutput() const; + void _stepOutput(const index_t it, const T resnorm, const T updatenorm) const; + + gsSparseMatrix m_sysmat; + + using Base::_computeForce; + using Base::_computeResidual; + using Base::_computeMass; + using Base::_computeMassInverse; + using Base::_computeDamping; + using Base::_computeJacobian; + +protected: + + using Base::m_solver; + + using Base::m_numIterations; + + using Base::m_options; + + +private: + template + typename std::enable_if<(_nonlinear==false), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; + + template + typename std::enable_if<(_nonlinear==true), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; +}; + +} // namespace gismo + +#ifndef GISMO_BUILD_LIB +#include GISMO_HPP_HEADER(gsDynamicImplicitEuler.hpp) +#endif diff --git a/src/gsDynamicSolvers/gsDynamicImplicitEuler.hpp b/src/gsDynamicSolvers/gsDynamicImplicitEuler.hpp new file mode 100644 index 0000000..c97ada0 --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicImplicitEuler.hpp @@ -0,0 +1,190 @@ +/** @file gsDynamicImplicitEuler.hpp + + @brief Performs the arc length method to solve a nonlinear equation system. + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-..., TU Delft) +*/ + +#pragma once + +#include +#include +#include + +namespace gismo +{ + +template +template +typename std::enable_if<(_nonlinear==false), gsStatus>::type +gsDynamicImplicitEuler::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + index_t N = U.rows(); + gsVector sol(2*N); + sol.topRows(N) = U; + sol.bottomRows(N) = V; + + typename gsBlockOp<>::Ptr Amat; + + Amat=gsBlockOp<>::make(2,2); + gsGMRes gmres(Amat); + + gsSparseMatrix eye(N,N); + eye.setIdentity(); + + gsVector F; + gsSparseMatrix M, C, K; + gsSparseMatrix Minv; + + // Computed at t=t0+dt + this->_computeMass(t+dt,M); + this->_computeForce(t+dt,F); + this->_computeDamping(U,t+dt,C); + this->_computeJacobian(U,t,K); + + // top-left + Amat->addOperator(0,0,gsIdentityOp::make(N) ); + // top-right + Amat->addOperator(0,1,makeMatrixOp( -dt*eye ) ); + // bottom-left + Amat->addOperator(1,0,makeMatrixOp( dt*K ) ); + // bottom-right + Amat->addOperator(1,1,makeMatrixOp( M + dt*C ) ); + + gsVector rhs(2*N); + rhs.topRows(N) = U; + rhs.bottomRows(N) = dt*F + M*V; + + gsMatrix tmpsol; + gmres.solve(rhs,tmpsol); + + U = tmpsol.topRows(N); + V = tmpsol.bottomRows(N); + this->_initOutput(); + this->_stepOutput(0,sol.norm(),0.); + + return gsStatus::Success; +} + + +template +template +typename std::enable_if<(_nonlinear==true), gsStatus>::type +gsDynamicImplicitEuler::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + index_t N = U.rows(); + gsVector sol(2*N); + sol.topRows(N) = U; + sol.bottomRows(N) = V; + + gsMatrix dsol; + + typename gsBlockOp<>::Ptr Amat; + + Amat=gsBlockOp<>::make(2,2); + gsGMRes gmres(Amat); + + gsSparseMatrix eye(N,N); + eye.setIdentity(); + + gsVector R; + gsSparseMatrix M, C, K; + + // Computed at t=t0+dt + this->_computeMass(t+dt,M); + this->_computeDamping(U,t+dt,C); + this->_computeResidual(U,t+dt,R); + + gsVector rhs(2*N); + rhs.topRows(N) = - dt * V; // same as below, but U-Uold=0 + rhs.bottomRows(N) = dt * C * V + dt*(-R); // same as below, but V-Vold=0 + + T tolU = m_options.getReal("TolU"); + T tolF = m_options.getReal("TolF"); + T updateNorm = 10.0*tolU; + T residualNorm = rhs.norm(); + T residualNorm0 = (residualNorm!=0) ? residualNorm : 1; + T solnorm, dsolnorm; + this->_initOutput(); + for (index_t numIterations = 0; numIterations < m_options.getInt("MaxIter"); ++numIterations) + { + // TODO: Quasi newton + if ((!m_options.getSwitch("Quasi")) || ((numIterations==0) || (numIterations % m_options.getInt("QuasiIterations") == 0))) + { + // Computed at t=t0+dt + this->_computeDamping(U,t+dt,C); + this->_computeJacobian(U,t+dt,K); + } + + Amat->addOperator(0,0,gsIdentityOp::make(N) ); + Amat->addOperator(0,1,makeMatrixOp(-dt*eye) ); + Amat->addOperator(1,0,makeMatrixOp(dt*K) ); + Amat->addOperator(1,1,makeMatrixOp(M + dt*C) ); + + rhs.topRows(N) = U - Uold - dt * V; + rhs.bottomRows(N) = M*(V - Vold) + dt * C * V + dt*(-R); + + gmres.solve(-rhs,dsol); + sol += dsol; + + solnorm = sol.norm(); + dsolnorm = dsol.norm(); + updateNorm = (solnorm != 0) ? dsolnorm / solnorm : dsolnorm; + + U = sol.topRows(N); + V = sol.bottomRows(N); + + this->_computeResidual(U,t,R); + residualNorm = rhs.norm() / residualNorm0; + + this->_stepOutput(numIterations,residualNorm,updateNorm); + + if ( (updateNorm +void gsDynamicImplicitEuler::_initOutput() const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< +void gsDynamicImplicitEuler::_stepOutput(const index_t it, const T resnorm, const T updatenorm) const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< + +#include +#include + +namespace gismo +{ + +/** + @brief Performs the arc length method to solve a nonlinear system of equations. + + \tparam T coefficient type + + \ingroup gsStructuralAnalysis +*/ +template +class gsDynamicNewmark : public gsDynamicBase +{ + typedef gsDynamicBase Base; + +protected: + + typedef typename gsStructuralAnalysisOps::Force_t Force_t; + typedef typename gsStructuralAnalysisOps::TForce_t TForce_t; + typedef typename gsStructuralAnalysisOps::Residual_t Residual_t; + typedef typename gsStructuralAnalysisOps::TResidual_t TResidual_t; + typedef typename gsStructuralAnalysisOps::Mass_t Mass_t; + typedef typename gsStructuralAnalysisOps::TMass_t TMass_t; + typedef typename gsStructuralAnalysisOps::Damping_t Damping_t; + typedef typename gsStructuralAnalysisOps::TDamping_t TDamping_t; + typedef typename gsStructuralAnalysisOps::Stiffness_t Stiffness_t; + typedef typename gsStructuralAnalysisOps::Jacobian_t Jacobian_t; + typedef typename gsStructuralAnalysisOps::TJacobian_t TJacobian_t; + +public: + + virtual ~gsDynamicNewmark() {}; + + /// Constructor + gsDynamicNewmark( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const Force_t & Force + ) + : + Base(Mass,Damping,Stiffness,Force) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicNewmark( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const TForce_t & TForce + ) + : + Base(Mass,Damping,Stiffness,TForce) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicNewmark( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const Residual_t & Residual + ) + : + Base(Mass,Damping,Jacobian,Residual) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicNewmark( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,Jacobian,TResidual) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicNewmark( + const Mass_t & Mass, + const Damping_t & Damping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,TJacobian,TResidual) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicNewmark( + const TMass_t & TMass, + const TDamping_t & TDamping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(TMass,TDamping,TJacobian,TResidual) + { + this->defaultOptions(); + } + + /// Set default options + virtual void defaultOptions() override; + +// General functions +protected: + + gsStatus _step(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const override + { + gsStatus status = gsStatus::NotStarted; + status = _step_impl<_NL>(t,dt,U,V,A); + return status; + } + + void _initOutput() const; + void _stepOutput(const index_t it, const T resnorm, const T updatenorm) const; + + gsSparseMatrix m_sysmat; + + using Base::_computeForce; + using Base::_computeResidual; + using Base::_computeMass; + using Base::_computeMassInverse; + using Base::_computeDamping; + using Base::_computeJacobian; + +protected: + + using Base::m_solver; + + using Base::m_numIterations; + + using Base::m_options; + + +private: + template + typename std::enable_if<(_nonlinear==false), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; + + template + typename std::enable_if<(_nonlinear==true), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; +}; + +} // namespace gismo + +#ifndef GISMO_BUILD_LIB +#include GISMO_HPP_HEADER(gsDynamicNewmark.hpp) +#endif diff --git a/src/gsDynamicSolvers/gsDynamicNewmark.hpp b/src/gsDynamicSolvers/gsDynamicNewmark.hpp new file mode 100644 index 0000000..c72895c --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicNewmark.hpp @@ -0,0 +1,176 @@ +/** @file gsDynamicNewmark.hpp + + @brief Performs the arc length method to solve a nonlinear equation system. + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-..., TU Delft) +*/ + +#pragma once + +namespace gismo +{ + +template +void gsDynamicNewmark::defaultOptions() +{ + Base::defaultOptions(); + m_options.addReal("alpha","Beta parameter for Newmark's method, such that 0 =< 2 beta =< 1",0.25); + m_options.addReal("delta","Beta parameter for Newmark's method, such that 0 =< gamma =< 1",0.50); +} + + +template +template +typename std::enable_if<(_nonlinear==false), gsStatus>::type +gsDynamicNewmark::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + gsVector F; + gsSparseMatrix M, C, K; + + T alpha = m_options.getReal("alpha"); + T delta = m_options.getReal("delta"); + + // Computed at t=t0+dt + this->_computeMass(t+dt,M); + this->_computeForce(t+dt,F); + + gsSparseMatrix lhs; + gsMatrix rhs; + // predictors + U = Uold + dt*Vold + Aold*(0.5 - alpha)*dt*dt; + V = Vold + Aold*(1 - delta)*dt; + + this->_computeDamping(U,t+dt,C); + this->_computeJacobian(U,t+dt,K); + + // lhs and rhs + lhs = M + delta*dt*C + dt*dt*alpha*K; + rhs = F - K*U - C * V; + + this->_initOutput(); + typename gsSparseSolver::uPtr solver; + solver = gsSparseSolver::get( m_options.getString("Solver") ); + solver->compute(lhs); + + A = solver->solve(rhs); + V += A*delta*dt; + U += A*alpha*dt*dt; + + this->_stepOutput(0,(F - K*U - M * A - C * V).norm(),U.norm()); + if (math::isinf(U.norm()) || math::isnan(U.norm())) + return gsStatus::NotConverged; + else + return gsStatus::Success; +} + + +template +template +typename std::enable_if<(_nonlinear==true), gsStatus>::type +gsDynamicNewmark::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + //https://ethz.ch/content/dam/ethz/special-interest/baug/ibk/structural-mechanics-dam/education/femII/presentation_05_dynamics_v3.pdf + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + gsVector R; + gsSparseMatrix M, C, K; + + T alpha = m_options.getReal("alpha"); + T delta = m_options.getReal("delta"); + + gsSparseMatrix lhs; + A.setZero(); + U = Uold + dt*Vold + Aold*(0.5 - alpha)*dt*dt + alpha*dt*dt*A; + V = Vold + Aold*(1 - delta)*dt + delta*dt*A; + // Computed at t=t0+dt + this->_computeMass(t+dt,M); + this->_computeDamping(U,t+dt,C); + this->_computeResidual(U,t+dt,R); + + gsMatrix rhs = R - C * V - M*(A); + + T tolU = m_options.getReal("TolU"); + T tolF = m_options.getReal("TolF"); + T updateNorm = 10.0*tolU; + T residualNorm = rhs.norm(); + T residualNorm0 = (residualNorm!=0) ? residualNorm : 1; + gsVector dA; + T Anorm, dAnorm; + this->_initOutput(); + typename gsSparseSolver::uPtr solver; + for (index_t numIterations = 0; numIterations < m_options.getInt("MaxIter"); ++numIterations) + { + if ((!m_options.getSwitch("Quasi")) || ((numIterations==0) || (numIterations % m_options.getInt("QuasiIterations") == 0))) + { + // Computed at t=t0+dt + this->_computeDamping(U,t+dt,C); + this->_computeJacobian(U,t+dt,K); + } + + lhs = M + delta*dt*C + dt*dt*alpha*K; + + solver = gsSparseSolver::get( m_options.getString("Solver") ); + solver->compute(lhs); + + dA = solver->solve(rhs); + A += dA; + V += dA*delta*dt; + U += dA*alpha*dt*dt; + + Anorm = A.norm(); + dAnorm = dA.norm(); + updateNorm = (Anorm!=0) ? dAnorm/Anorm : dAnorm; + + this->_computeResidual(U,t+dt,R); + rhs = R - C*V - M*A; + residualNorm = rhs.norm() / residualNorm0; + + this->_stepOutput(numIterations,residualNorm,updateNorm); + + if ( (updateNorm +void gsDynamicNewmark::_initOutput() const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< +void gsDynamicNewmark::_stepOutput(const index_t it, const T resnorm, const T updatenorm) const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< + +#include + +namespace gismo +{ + +/** + @brief Performs the arc length method to solve a nonlinear system of equations. + + \tparam T coefficient type + + \ingroup gsStructuralAnalysis +*/ +template +class gsDynamicRK4 : public gsDynamicBase +{ + typedef gsDynamicBase Base; + +protected: + + typedef typename gsStructuralAnalysisOps::Force_t Force_t; + typedef typename gsStructuralAnalysisOps::TForce_t TForce_t; + typedef typename gsStructuralAnalysisOps::Residual_t Residual_t; + typedef typename gsStructuralAnalysisOps::TResidual_t TResidual_t; + typedef typename gsStructuralAnalysisOps::Mass_t Mass_t; + typedef typename gsStructuralAnalysisOps::TMass_t TMass_t; + typedef typename gsStructuralAnalysisOps::Damping_t Damping_t; + typedef typename gsStructuralAnalysisOps::TDamping_t TDamping_t; + typedef typename gsStructuralAnalysisOps::Stiffness_t Stiffness_t; + typedef typename gsStructuralAnalysisOps::Jacobian_t Jacobian_t; + typedef typename gsStructuralAnalysisOps::TJacobian_t TJacobian_t; + +public: + + virtual ~gsDynamicRK4() {}; + + /// Constructor + gsDynamicRK4( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const Force_t & Force + ) + : + Base(Mass,Damping,Stiffness,Force) + {} + + /// Constructor + gsDynamicRK4( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const TForce_t & TForce + ) + : + Base(Mass,Damping,Stiffness,TForce) + {} + + /// Constructor + gsDynamicRK4( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const Residual_t & Residual + ) + : + Base(Mass,Damping,Jacobian,Residual) + {} + + /// Constructor + gsDynamicRK4( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,Jacobian,TResidual) + {} + + /// Constructor + gsDynamicRK4( + const Mass_t & Mass, + const Damping_t & Damping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,TJacobian,TResidual) + {} + + /// Constructor + gsDynamicRK4( + const TMass_t & TMass, + const TDamping_t & TDamping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(TMass,TDamping,TJacobian,TResidual) + {} + +// General functions +protected: + + gsStatus _step(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const + { + gsStatus status = gsStatus::NotStarted; + status = _step_impl<_NL>(t,dt,U,V,A); + return status; + } + + void _initOutput() const; + void _stepOutput(const index_t it, const T resnorm, const T updatenorm) const; + + gsSparseMatrix m_sysmat; + + using Base::_computeForce; + using Base::_computeResidual; + using Base::_computeMass; + using Base::_computeMassInverse; + using Base::_computeDamping; + using Base::_computeJacobian; + +protected: + + using Base::m_solver; + + using Base::m_numIterations; + + using Base::m_options; + + +private: + template + typename std::enable_if<(_nonlinear==false), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; + + template + typename std::enable_if<(_nonlinear==true), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; +}; + +} // namespace gismo + +#ifndef GISMO_BUILD_LIB +#include GISMO_HPP_HEADER(gsDynamicRK4.hpp) +#endif diff --git a/src/gsDynamicSolvers/gsDynamicRK4.hpp b/src/gsDynamicSolvers/gsDynamicRK4.hpp new file mode 100644 index 0000000..d440e66 --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicRK4.hpp @@ -0,0 +1,183 @@ +/** @file gsDynamicRK4.hpp + + @brief Performs the arc length method to solve a nonlinear equation system. + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-..., TU Delft) +*/ + +#pragma once + +namespace gismo +{ +template +template +typename std::enable_if<(_nonlinear==false), gsStatus>::type +gsDynamicRK4::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + /* + */ + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + index_t N = U.rows(); + gsVector sol(2*N); + sol.topRows(N) = Uold; + sol.bottomRows(N) = Vold; + + gsVector F, R; //R is the residual vector, and the new _computeForce give the inline factor + gsSparseMatrix M, Minv, C, K; + + // Computed at t=t0 + this->_computeMass(t,M); + this->_computeMassInverse(M,Minv); + // this->_computeForce(t,F); + this->_computeDamping(U,t,C); //C is damping + this->_computeJacobian(U,t,K); + + // this->_initOutput(); + // Initialize parameters for RK4 + gsVector k1(2*N), k2(2*N), k3(2*N), k4(2*N); + gsVector Utmp, Vtmp; + + //Step1 (calculate k1) + _computeForce(t, F); + R = F - K * Uold; + k1.topRows(N) = Vold; + k1.bottomRows(N) = Minv * (R - C * Vold); + + //Step2 (calculate k2) + Utmp = Uold + dt/2. * k1.topRows(N); + Vtmp = Vold + dt/2. * k1.bottomRows(N); + _computeForce(t + dt/2.,F); + R = F - K * Utmp; + k2.topRows(N) = Vtmp; + k2.bottomRows(N) = Minv * ( R - C * Vtmp); + + //Step3 (calculate k3) + Utmp = Uold + dt/2. * k2.topRows(N); + Vtmp = Vold + dt/2. * k2.bottomRows(N); + _computeForce(t + dt/2., F); + R = F - K * Utmp; + k3.topRows(N) = Vtmp; + k3.bottomRows(N) = Minv * ( R - C * Vtmp); + + //Step4 (calculate k4) + Utmp = Uold + dt/2. * k3.topRows(N); + Vtmp = Vold + dt/2. * k3.bottomRows(N); + _computeForce(t + dt/2., F); + R = F - K * Utmp; + k4.topRows(N) = Vtmp; + k4.bottomRows(N) = Minv * ( R - C * Vtmp); + + sol += 1./6 * dt * (k1 + 2.*k2 + 2.*k3 + k4); + + U = std::move(sol.topRows(N)); //Use std move to update vectors to avoid unnecessary copying + V = std::move(sol.bottomRows(N)); + + if (math::isinf(sol.norm()) || math::isnan(sol.norm())) + return gsStatus::NotConverged; + else + return gsStatus::Success; +} + + +template +template +typename std::enable_if<(_nonlinear==true), gsStatus>::type +gsDynamicRK4::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + /* + */ + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + index_t N = U.rows(); + gsVector sol(2*N); + sol.topRows(N) = Uold; + sol.bottomRows(N) = Vold; + + gsVector F, R; //R is the residual vector, and the new _computeForce give the inline factor + gsSparseMatrix M, Minv, C, K; + + // Computed at t=t0 + this->_computeMass(t,M); + this->_computeMassInverse(M,Minv); + // this->_computeForce(t,F); + this->_computeDamping(U,t,C); //C is damping + this->_computeJacobian(U,t,K); + + // this->_initOutput(); + // Initialize parameters for RK4 + gsVector k1(2*N), k2(2*N), k3(2*N), k4(2*N); + gsVector Utmp, Vtmp; + + //Step1 (calculate k1) + _computeResidual(Uold, t, R); + k1.topRows(N) = Vold; + k1.bottomRows(N) = Minv * (R - C * Vold); + + //Step2 (calculate k2) + Utmp = Uold + dt/2. * k1.topRows(N); + Vtmp = Vold + dt/2. * k1.bottomRows(N); + _computeResidual(Utmp,t + dt/2.,R); + k2.topRows(N) = Vtmp; + k2.bottomRows(N) = Minv * ( R - C * Vtmp); + + //Step3 (calculate k3) + Utmp = Uold + dt/2. * k2.topRows(N); + Vtmp = Vold + dt/2. * k2.bottomRows(N); + _computeResidual(Utmp,t + dt/2.,R); + k3.topRows(N) = Vtmp; + k3.bottomRows(N) = Minv * ( R - C * Vtmp); + + //Step4 (calculate k4) + Utmp = Uold + dt/2. * k3.topRows(N); + Vtmp = Vold + dt/2. * k3.bottomRows(N); + _computeResidual(Utmp,t + dt/2.,R); + k4.topRows(N) = Vtmp; + k4.bottomRows(N) = Minv * ( R - C * Vtmp); + + sol += 1./6 * dt * (k1 + 2.*k2 + 2.*k3 + k4); + + U = std::move(sol.topRows(N)); //Use std move to update vectors to avoid unnecessary copying + V = std::move(sol.bottomRows(N)); + + if (math::isinf(sol.norm()) || math::isnan(sol.norm())) + return gsStatus::NotConverged; + else + return gsStatus::Success; +} + +template +void gsDynamicRK4::_initOutput() const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< +void gsDynamicRK4::_stepOutput(const index_t it, const T resnorm, const T updatenorm) const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< + +#include +#include +#include + +namespace gismo +{ + +/** + @brief Performs the arc length method to solve a nonlinear system of equations. + + \tparam T coefficient type + + \ingroup gsStructuralAnalysis +*/ +template +class gsDynamicWilson : public gsDynamicNewmark +{ + typedef gsDynamicNewmark Base; + +protected: + + typedef typename gsStructuralAnalysisOps::Force_t Force_t; + typedef typename gsStructuralAnalysisOps::TForce_t TForce_t; + typedef typename gsStructuralAnalysisOps::Residual_t Residual_t; + typedef typename gsStructuralAnalysisOps::TResidual_t TResidual_t; + typedef typename gsStructuralAnalysisOps::Mass_t Mass_t; + typedef typename gsStructuralAnalysisOps::TMass_t TMass_t; + typedef typename gsStructuralAnalysisOps::Damping_t Damping_t; + typedef typename gsStructuralAnalysisOps::TDamping_t TDamping_t; + typedef typename gsStructuralAnalysisOps::Stiffness_t Stiffness_t; + typedef typename gsStructuralAnalysisOps::Jacobian_t Jacobian_t; + typedef typename gsStructuralAnalysisOps::TJacobian_t TJacobian_t; + +public: + + virtual ~gsDynamicWilson() {}; + + /// Constructor + gsDynamicWilson( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const Force_t & Force + ) + : + Base(Mass,Damping,Stiffness,Force) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicWilson( + const Mass_t & Mass, + const Damping_t & Damping, + const Stiffness_t & Stiffness, + const TForce_t & TForce + ) + : + Base(Mass,Damping,Stiffness,TForce) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicWilson( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const Residual_t & Residual + ) + : + Base(Mass,Damping,Jacobian,Residual) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicWilson( + const Mass_t & Mass, + const Damping_t & Damping, + const Jacobian_t & Jacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,Jacobian,TResidual) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicWilson( + const Mass_t & Mass, + const Damping_t & Damping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(Mass,Damping,TJacobian,TResidual) + { + this->defaultOptions(); + } + + /// Constructor + gsDynamicWilson( + const TMass_t & TMass, + const TDamping_t & TDamping, + const TJacobian_t & TJacobian, + const TResidual_t & TResidual + ) + : + Base(TMass,TDamping,TJacobian,TResidual) + { + this->defaultOptions(); + } + + /// Set default options + virtual void defaultOptions() override; + +// General functions +protected: + + gsStatus _step(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const override + { + gsStatus status = gsStatus::NotStarted; + status = _step_impl<_NL>(t,dt,U,V,A); + return status; + } + + void _initOutput() const; + void _stageOutput(index_t stage) const; + void _stepOutput(const index_t it, const T resnorm, const T updatenorm) const; + + gsSparseMatrix m_sysmat; + + using Base::_computeForce; + using Base::_computeResidual; + using Base::_computeMass; + using Base::_computeMassInverse; + using Base::_computeDamping; + using Base::_computeJacobian; + +protected: + + using Base::m_solver; + + using Base::m_numIterations; + + using Base::m_options; + + +private: + template + typename std::enable_if<(_nonlinear==false), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; + + template + typename std::enable_if<(_nonlinear==true), gsStatus>::type + _step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const; +}; + +} // namespace gismo + +#ifndef GISMO_BUILD_LIB +#include GISMO_HPP_HEADER(gsDynamicWilson.hpp) +#endif diff --git a/src/gsDynamicSolvers/gsDynamicWilson.hpp b/src/gsDynamicSolvers/gsDynamicWilson.hpp new file mode 100644 index 0000000..40c3886 --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicWilson.hpp @@ -0,0 +1,120 @@ +/** @file gsDynamicWilson.hpp + + @brief Performs the arc length method to solve a nonlinear equation system. + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst (2019-..., TU Delft) +*/ + +#pragma once + +namespace gismo +{ + +template +void gsDynamicWilson::defaultOptions() +{ + Base::defaultOptions(); + m_options.addReal("alpha","Beta parameter for Wilson's method, such that 0 =< 2 beta =< 1",0.25); + m_options.addReal("delta","Beta parameter for Wilson's method, such that 0 =< delta =< 1",0.50); + m_options.addReal("gamma","Beta parameter for Wilson's method, such that 0 =< gamma =< 1",1.4); +} + + +template +template +typename std::enable_if<(_nonlinear==false), gsStatus>::type +gsDynamicWilson::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + T alpha = m_options.getReal("alpha"); + T delta = m_options.getReal("delta"); + T gamma = m_options.getReal("gamma"); + + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + gsVector Utmp = U; + gsVector Vtmp = V; + gsVector Atmp = A; + + // Perform a Newmark step with DT = gamma*dt + Base::_step(t,gamma*dt,Utmp,Vtmp,Atmp); + + // Compute solutions + A = (1-1/gamma)*Aold + 1/gamma*Atmp; + V += A*delta*dt; + U += A*alpha*dt*dt; + + return gsStatus::Success; +} + + +template +template +typename std::enable_if<(_nonlinear==true), gsStatus>::type +gsDynamicWilson::_step_impl(const T t, const T dt, gsVector & U, gsVector & V, gsVector & A) const +{ + T alpha = m_options.getReal("alpha"); + T delta = m_options.getReal("delta"); + T gamma = m_options.getReal("gamma"); + + gsVector Uold = U; + gsVector Vold = V; + gsVector Aold = A; + + gsVector Utmp = U; + gsVector Vtmp = V; + gsVector Atmp = A; + + /// Perform a Newmark step with DT = gamma*dt + Base::_step(t,gamma*dt,Utmp,Vtmp,Atmp); + + // Compute solutions (see Wilson et al. (1973) NONLINEAR DYNAMIC ANALYSIS OF COMPLEX STRUCTURES) + A = (1-1/gamma)*Aold + 1/gamma*Atmp; + U = Uold + dt*Vold + Aold*(0.5 - alpha)*dt*dt + alpha*dt*dt*A; + V = Vold + Aold*(1 - delta)*dt + delta*dt*A; + + return gsStatus::Success; +} + +template +void gsDynamicWilson::_stageOutput(index_t stage) const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"Stage "< +void gsDynamicWilson::_initOutput() const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< +void gsDynamicWilson::_stepOutput(const index_t it, const T resnorm, const T updatenorm) const +{ + if (m_options.getSwitch("Verbose")) + { + gsInfo<<"\t"; + gsInfo< +#include +#include + +namespace gismo +{ + +/** + @brief Performs the arc length method to solve a nonlinear system of equations. + + \tparam T coefficient type + + \ingroup gsStructuralAnalysis +*/ +template +class gsDynamicXBraid : public gsXBraid< gsMatrix > +{ + +public: + typedef typename std::function&,const gsVector&,const gsVector&)> callback_type; + + virtual ~gsDynamicXBraid() {}; + + /// Constructor + gsDynamicXBraid( + const gsDynamicBase * solver, + const gsMpiComm& comm, + const T& tstart, + const T& tend, + const index_t numDofs, + index_t numSteps = 10) + : + gsXBraid< gsMatrix >::gsXBraid(comm, tstart, tend, (int)numSteps), + m_solver(solver), + m_numDofs(numDofs) + { + defaultOptions(); + m_callback = [](const index_t&,const T&,const gsVector&,const gsVector&,const gsVector&){return;}; // set empty callback + } + +// Member functions +public: + void defaultOptions() + { + // see gsXBraid/filedata/pde/heat2d_square_ibvp1.xml + m_options.addInt("CFactor", "Coarsening factor of the parallel-in-time multigrid solver", 2); + m_options.addInt("access", "Access level (never [=0], =after finished [=1(default)], each iteration [=2]", 1); + m_options.addInt("maxIter", "Maximum iteration numbers of the parallel-in-time multigrid solver", 100); + m_options.addInt("maxLevel", "Maximum numbers of parallel-in-time multigrid levels", 30); + m_options.addInt("minCLevel", "Minimum level of the parallel-in-time multigrid solver", 2); + m_options.addInt("norm", "Temporal norm of the parallel-in-time multigrid solver (1-norm [=1], 2-norm [=2(default)], inf-norm [=3])", 2); + m_options.addInt("numFMG", "Number of full multigrid steps of the parallel-in-time multigrid solver", 1); + m_options.addInt("numFMGVcyc", "Number of full multigrid V-cycles of the parallel-in-time multigrid solver", 1); + m_options.addInt("numMaxRef", "Maximum number of refinements of the parallel-in-time multigrid solver", 1); + m_options.addInt("numRelax", "Number of relaxation steps of the parallel-in-time multigrid solver", 1); + m_options.addInt("numStorage", "Number of storage of the parallel-in-time multigrid solver", -1); + m_options.addInt("print", "Print level (no output [=0], runtime inforation [=1], run statistics [=2(default)], debug [=3])", 2); + m_options.addReal("absTol", "Absolute tolerance of the parallel-in-time multigrid solver", 1e-6); + m_options.addReal("relTol", "Relative tolerance of the parallel-in-time multigrid solver", 1e-3); + m_options.addSwitch("fmg", "Perform full multigrid (default is off)", 0); + m_options.addSwitch("incrMaxLevels", "Increase the maximum number of parallel-in-time multigrid levels after performing a refinement (default is off)", 0); + m_options.addSwitch("periodic", "Periodic time grid (default is off)", 0); + m_options.addSwitch("refine", "Perform refinement in time (default off)", 0); + m_options.addSwitch("sequential", "Set the initial guess of the parallel-in-time multigrid solver as the sequential time stepping solution (default is off)", 0); + m_options.addSwitch("skip", "Skip all work on the first down cycle of the parallel-in-time multigrid solver (default on)", 1); + m_options.addSwitch("spatial", "Perform spatial coarsening and refinement (default is off)", 0); + m_options.addSwitch("tol", "Tolerance type (absolute [=true], relative [=false(default)]", 0); + + + m_options.addSwitch("extraVerbose", "Extra verbosity", 0); + } + + void initialize() + { + // see gsXBraid/examples/xbraid_heatEquation_example.cpp + this->SetCFactor(m_options.getInt("CFactor")); + this->SetMaxIter(m_options.getInt("maxIter")); + this->SetMaxLevels(m_options.getInt("maxLevel")); + this->SetMaxRefinements(m_options.getInt("numMaxRef")); + this->SetMinCoarse(m_options.getInt("minCLevel")); + this->SetNFMG(m_options.getInt("numFMG")); + this->SetNFMGVcyc(m_options.getInt("numFMGVcyc")); + this->SetNRelax(m_options.getInt("numRelax")); + this->SetAccessLevel(m_options.getInt("access")); + this->SetPrintLevel(m_options.getInt("print")); + this->SetStorage(m_options.getInt("numStorage")); + this->SetTemporalNorm(m_options.getInt("norm")); + + if (m_options.getSwitch("fmg")) this->SetFMG(); + if (m_options.getSwitch("incrMaxLevels")) this->SetIncrMaxLevels(); + if (m_options.getSwitch("periodic")) this->SetPeriodic(1); else this->SetPeriodic(0); + if (m_options.getSwitch("refine")) this->SetRefine(1); else this->SetRefine(0); + if (m_options.getSwitch("sequential")) this->SetSeqSoln(1); else this->SetSeqSoln(0); + if (m_options.getSwitch("skip")) this->SetSkip(1); else this->SetSkip(0); + if (m_options.getSwitch("spatial")) this->SetSpatialCoarsenAndRefine(); + if (m_options.getSwitch("tol")) this->SetAbsTol(m_options.getReal("absTol")); + else this->SetRelTol(m_options.getReal("relTol")); + } + + /// Initializes a vector + braid_Int Init(braid_Real t, braid_Vector *u_ptr) override + { + gsMatrix* u = new gsMatrix(3*m_numDofs, 1); + + // Does this mean zero displacements? + u->setZero(); + + if (m_solver->solutionU().rows()==m_numDofs) + u->col(0).segment(0 ,m_numDofs) = m_solver->solutionU(); + if (m_solver->solutionV().rows()==m_numDofs) + u->col(0).segment(m_numDofs ,m_numDofs) = m_solver->solutionV(); + if (m_solver->solutionA().rows()==m_numDofs) + u->col(0).segment(2*m_numDofs,m_numDofs) = m_solver->solutionA(); + + *u_ptr = (braid_Vector) u; + return braid_Int(0); + } + + gsOptionList & options() { return m_options; } + void setOptions(gsOptionList options) {m_options.update(options,gsOptionList::addIfUnknown); }; + + /// Performs a single step of the parallel-in-time multigrid + braid_Int Step(braid_Vector u, braid_Vector ustop, braid_Vector fstop, BraidStepStatus &status) override + { + gsVector* u_ptr = (gsVector*) u; + // gsMatrix* ustop_ptr = (gsMatrix*) ustop; // the guess is not used + + // XBraid forcing + if (fstop != NULL) + { + gsVector* fstop_ptr = (gsVector*) fstop; + *u_ptr += *fstop_ptr; + } + + // Get time step information + std::pair time = static_cast(status).timeInterval(); + if (m_options.getSwitch("extraVerbose")) gsInfo<<"Solving interval ["<(status).level()<<")\n"; + T t = time.first; + T dt = time.second - time.first; + + // Solve time step + gsVector U = (*u_ptr).segment(0 ,m_numDofs); + gsVector V = (*u_ptr).segment(m_numDofs ,m_numDofs); + gsVector A = (*u_ptr).segment(2*m_numDofs,m_numDofs); + + gsStatus stepStatus = m_solver->step(t,dt,U,V,A); + + u_ptr->segment(0 ,m_numDofs) = U; + u_ptr->segment(m_numDofs ,m_numDofs) = V; + u_ptr->segment(2*m_numDofs,m_numDofs) = A; + + // Carry out adaptive refinement in time + if (static_cast(status).level() == 0) + { + if (stepStatus==gsStatus::Success) + { + braid_Real error = static_cast(status).error(); + if (error != braid_Real(-1.0)) + { + braid_Int rfactor = (braid_Int) std::ceil( std::sqrt( error / 1e-3) ); + status.SetRFactor(rfactor); + } + else + status.SetRFactor(1); + } + // Refine if solution interval failed + else + { + if (m_options.getSwitch("extraVerbose")) gsInfo<<"Step "<<(static_cast(status)).timeIndex()<<" did not converge"; + status.SetRFactor((braid_Int)2); + } + } + return braid_Int(0); + } + + /// Computes the spatial norm of the given vector + braid_Int SpatialNorm( braid_Vector u, + braid_Real *norm_ptr) override + { + gsVector* u_ptr = (gsVector*) u; + *norm_ptr = u_ptr->segment(0,m_numDofs).norm(); // Displacement-based norm + // *norm_ptr = u_ptr->norm(); + return braid_Int(0); + } + + /// Sets the size of the MPI communication buffer + braid_Int BufSize(braid_Int *size_ptr, BraidBufferStatus &status) override + { + *size_ptr = sizeof(T)*(m_numDofs*3+2); // +2 comes from rows, cols of the solution vector u. + return braid_Int(0); + } + + void setCallback(callback_type callback) const {m_callback = callback;} + + /// Handles access for input/output + braid_Int Access(braid_Vector u, BraidAccessStatus &status) override + { + gsVector* u_ptr = (gsVector*) u; + m_callback((index_t) static_cast(status).timeIndex(), + (T) static_cast(status).time(), + (gsVector)(*u_ptr).segment(0 ,m_numDofs), + (gsVector)(*u_ptr).segment(m_numDofs ,m_numDofs), + (gsVector)(*u_ptr).segment(2*m_numDofs,m_numDofs) + ); + return braid_Int(0); + } + + /// Performs spatial coarsening + /* + NOTE: This routine is not implemented. How to do it: + 1. Make the Coarsen and Refine routines virtual + 2. Within the example, overload this class, and define Coarsen and Refine for the problem + 3. Update the m_numDoFs within!! + */ + braid_Int Coarsen(braid_Vector fu, braid_Vector *cu_ptr, BraidCoarsenRefStatus &status) override + { + gsMatrix *fu_ptr = (gsMatrix*) fu; + gsMatrix* cu = new gsMatrix(); + *cu = *fu_ptr; + *cu_ptr = (braid_Vector) cu; + return braid_Int(0); + } + + // Performs spatial refinement + braid_Int Refine(braid_Vector cu, braid_Vector *fu_ptr, BraidCoarsenRefStatus &status) override + { + gsMatrix *cu_ptr = (gsMatrix*) cu; + gsMatrix* fu = new gsMatrix(); + *fu = *cu_ptr; + *fu_ptr = (braid_Vector) fu; + return braid_Int(0); + } + + +// Class members +protected: + + const gsDynamicBase * m_solver; + index_t m_numDofs; + gsOptionList m_options; + mutable callback_type m_callback; + +}; + +} // namespace gismo + +#ifndef GISMO_BUILD_LIB +#include GISMO_HPP_HEADER(gsDynamicImplicitEuler.hpp) +#endif diff --git a/src/gsDynamicSolvers/gsDynamicXBraid_.cpp b/src/gsDynamicSolvers/gsDynamicXBraid_.cpp new file mode 100644 index 0000000..34ac986 --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamicXBraid_.cpp @@ -0,0 +1,12 @@ +#include + +#ifdef gsXbraid_ENABLED +#include +#endif + +namespace gismo +{ +#ifdef gsXbraid_ENABLED + CLASS_TEMPLATE_INST gsDynamicXBraid; +#endif +} diff --git a/src/gsDynamicSolvers/gsDynamic_.cpp b/src/gsDynamicSolvers/gsDynamic_.cpp new file mode 100644 index 0000000..d5b28ea --- /dev/null +++ b/src/gsDynamicSolvers/gsDynamic_.cpp @@ -0,0 +1,56 @@ +#include + +#include +#include + +#include +#include + +#include +#include + +#include +#include + +// #include +// #include + +#include +#include + +#include +#include + +#include +#include + +namespace gismo +{ + + enum struct gsStatus; + template struct gsStructuralAnalysisOps; + + CLASS_TEMPLATE_INST gsDynamicBase; + + CLASS_TEMPLATE_INST gsDynamicExplicitEuler; + CLASS_TEMPLATE_INST gsDynamicExplicitEuler; + + CLASS_TEMPLATE_INST gsDynamicImplicitEuler; + CLASS_TEMPLATE_INST gsDynamicImplicitEuler; + + CLASS_TEMPLATE_INST gsDynamicRK4; + CLASS_TEMPLATE_INST gsDynamicRK4; + + // CLASS_TEMPLATE_INST gsDynamicCentralDifference; + // CLASS_TEMPLATE_INST gsDynamicCentralDifference; + + CLASS_TEMPLATE_INST gsDynamicNewmark; + CLASS_TEMPLATE_INST gsDynamicNewmark; + + CLASS_TEMPLATE_INST gsDynamicWilson; + CLASS_TEMPLATE_INST gsDynamicWilson; + + CLASS_TEMPLATE_INST gsDynamicBathe; + CLASS_TEMPLATE_INST gsDynamicBathe; + +} diff --git a/src/gsDynamicSolvers/gsTimeIntegrator.h b/src/gsDynamicSolvers/gsTimeIntegrator.h deleted file mode 100644 index afd7c26..0000000 --- a/src/gsDynamicSolvers/gsTimeIntegrator.h +++ /dev/null @@ -1,493 +0,0 @@ -/** @file gsTimeIntegrator.h - - @brief Provides temporal solvers for structural analysis problems - - This file is part of the G+Smo library. - - This Source Code Form is subject to the terms of the Mozilla Public - License, v. 2.0. If a copy of the MPL was not distributed with this - file, You can obtain one at http://mozilla.org/MPL/2.0/. - - Author(s): H.M. Verhelst (2019-..., TU Delft) -*/ - - -#pragma once - -#include -#include - -namespace gismo -{ -/** - @brief Provides temporal solvers for structural analysis problems - - \tparam T coefficient type - - \ingroup gsTimeIntegrator -*/ -template -class gsTimeIntegrator -{ - - typedef typename gsStructuralAnalysisOps::TForce_t TForce_t; - typedef typename gsStructuralAnalysisOps::TResidual_t TResidual_t; - typedef typename gsStructuralAnalysisOps::Jacobian_t Jacobian_t; - typedef typename gsStructuralAnalysisOps::dJacobian_t dJacobian_t; - -public: - - - /** - * @brief Constructor - * - * @param[in] Mass The mass matrix - * @param[in] Damp The damping matrix - * @param[in] Stif The linear stiffness matrix - * @param[in] Force The time-dependent force - * @param[in] dt The time step - */ - gsTimeIntegrator( const gsSparseMatrix &Mass, - const gsSparseMatrix &Damp, - const gsSparseMatrix &Stif, - const TForce_t &Force, - const T dt - ) - : m_mass(Mass), - m_damp(Damp), - m_stif(Stif), - m_jacobian(nullptr), - m_djacobian(nullptr), - m_forceFun(Force), - m_dt(dt) - { - m_dofs = m_mass.cols(); - m_method = "ExplEuler"; - m_NL = false; - this->initializeCoefficients(); - this->initiate(); - } - - /** - * @brief Constructor - * - * @param[in] Mass The mass matrix - * @param[in] Stif The linear stiffness matrix - * @param[in] Force The time-dependent force - * @param[in] dt The time step - */ - gsTimeIntegrator( const gsSparseMatrix &Mass, - const gsSparseMatrix &Stif, - const TForce_t &Force, - const T dt - ) - : m_mass(Mass), - m_stif(Stif), - m_jacobian(nullptr), - m_djacobian(nullptr), - m_forceFun(Force), - m_dt(dt) - { - m_dofs = m_mass.cols(); - m_damp = gsSparseMatrix(m_dofs,m_dofs); - m_method = "ExplEuler"; - m_NL = false; - this->initializeCoefficients(); - this->initiate(); - m_solver = gsSparseSolver::get( "LU" ); // todo: make proper options - } - -private: - - /** - * @brief Constructor - * - * @param[in] Mass The mass matrix - * @param[in] Damp The damping matrix - * @param[in] Residual The residual vector - * @param[in] dt The time step - */ - gsTimeIntegrator( const gsSparseMatrix &Mass, - const gsSparseMatrix &Damp, - const TResidual_t &Residual, - const T dt - ) - : m_mass(Mass), - m_damp(Damp), - m_jacobian(nullptr), - m_residualFun(Residual), - m_dt(dt) - { - m_dofs = m_mass.cols(); - m_method = "Newmark"; - m_NL = true; - this->initializeCoefficients(); - this->initiate(); - m_solver = gsSparseSolver::get( "LU" ); // todo: make proper options - } - - /** - * @brief Constructor - * - * @param[in] Mass The mass matrix - * @param[in] Residual The residual vector - * @param[in] dt The time step - */ - gsTimeIntegrator( const gsSparseMatrix &Mass, - const TResidual_t &Residual, - const T dt - ) - : m_mass(Mass), - m_jacobian(nullptr), - m_residualFun(Residual), - m_dt(dt) - { - m_dofs = m_mass.cols(); - m_method = "Newmark"; - m_NL = true; - m_damp = gsSparseMatrix(m_dofs,m_dofs); - this->initializeCoefficients(); - this->initiate(); - m_solver = gsSparseSolver::get( "LU" ); // todo: make proper options - } - - -public: - - /** - * @brief Constructor - * - * @param[in] Mass The mass matrix - * @param[in] Damp The damping matrix - * @param[in] Jacobian The Jacobian matrix - * @param[in] Residual The residual vector - * @param[in] dt The time step - */ - gsTimeIntegrator( const gsSparseMatrix &Mass, - const gsSparseMatrix &Damp, - const Jacobian_t &Jacobian, - const TResidual_t &Residual, - const T dt - ) - : gsTimeIntegrator(Mass,Damp,Residual,dt) - { - m_jacobian = Jacobian; - m_djacobian = [this](gsVector const & x, gsVector const & dx, gsSparseMatrix & m) -> bool - { - return m_jacobian(x,m); - }; - m_solver = gsSparseSolver::get( "LU" ); // todo: make proper options - } - - - /** - * @brief Constructor - * - * @param[in] Mass The mass matrix - * @param[in] Jacobian The Jacobian matrix - * @param[in] Residual The residual vector - * @param[in] dt The time step - */ - gsTimeIntegrator( const gsSparseMatrix &Mass, - const Jacobian_t &Jacobian, - const TResidual_t &Residual, - const T dt - ) - : gsTimeIntegrator(Mass,Residual,dt) - { - m_jacobian = Jacobian; - m_djacobian = [this](gsVector const & x, gsVector const & dx, gsSparseMatrix & m) -> bool - { - return m_jacobian(x,m); - }; - m_solver = gsSparseSolver::get( "LU" ); // todo: make proper options - } - - - /** - * @brief Constructor - * - * @param[in] Mass The mass matrix - * @param[in] Damp The damping matrix - * @param[in] dJacobian The Jacobian matrix taking the solution as first argument and the update as second - * @param[in] Residual The residual vector - * @param[in] dt The time step - */ - gsTimeIntegrator( const gsSparseMatrix &Mass, - const gsSparseMatrix &Damp, - const dJacobian_t &dJacobian, - const TResidual_t &Residual, - const T dt - ) - : gsTimeIntegrator(Mass,Damp,Residual,dt) - { - m_djacobian = dJacobian; - m_solver = gsSparseSolver::get( "LU" ); // todo: make proper options - } - - - /** - * @brief Constructor - * - * @param[in] Mass The mass matrix - * @param[in] dJacobian The Jacobian matrix taking the solution as first argument and the update as second - * @param[in] Residual The residual vector - * @param[in] dt The time step - */ - gsTimeIntegrator( const gsSparseMatrix &Mass, - const dJacobian_t &dJacobian, - const TResidual_t &Residual, - const T dt - ) - : gsTimeIntegrator(Mass,Residual,dt) - { - m_djacobian = dJacobian; - m_damp = gsSparseMatrix(m_dofs,m_dofs); - m_solver = gsSparseSolver::get( "LU" ); // todo: make proper options - } - - -public: - /// Performs a step - gsStatus step(); - - /// Constructs the system - void constructSystem(); - -public: - - /// Tells if the Newton method converged - bool converged() const {return m_converged;} - - /// Returns the number of Newton iterations performed - index_t numIterations() const { return m_numIterations;} - - /// Returns the tolerance value used - T tolerance() const {return m_tolerance;} - - /// Returns the error after solving the nonlinear system - T residue() const {return m_residue;} - - /// Set the maximum number of Newton iterations allowed - void setMaxIterations(index_t nIter) {m_maxIterations = nIter;} - - /// Set the tolerance for convergence - void setTolerance(T tol) {m_tolerance = tol;} - - /// Sets the time - void setTime(T time) {m_t = time; } - - /// Sets the time - void setTimeStep(T dt) {m_dt = dt; } - /// Gets the time - T currentTime() const {return m_t; } - - /// Set mass matrix - void setMassMatrix(const gsSparseMatrix& Mass) - { - m_mass = Mass; - m_dofs = m_mass.cols(); - } - /// Set damping matrix - void setDampingMatrix(const gsSparseMatrix& Damp) - { - m_damp = Damp; - m_dofs = m_damp.cols(); - } - /// Reset damping matrix - void resetDampingMatrix() - { - m_damp = gsSparseMatrix(m_dofs,m_dofs); - } - /// Set stiffness matrix - void setStiffnessMatrix(const gsSparseMatrix& Stif) - { - m_stif = Stif; - m_dofs = m_mass.cols(); - } - /// Set Jacobian matrix - void setJacobian(Jacobian_t &Jacobian) - { - m_jacobian = Jacobian; - m_djacobian = [this](gsVector const & x, gsVector const & dx, gsSparseMatrix & m) - { - return m_jacobian(x,m); - }; - m_dofs = m_mass.cols(); - } - /// Set the Jacobian matrix taking the solution in the first argument and the update as second - void setJacobian(dJacobian_t &dJacobian) - { - m_djacobian = dJacobian; - } - /// Resets the solution - void resetSolution(); - - /// Set the displacement - void setDisplacement(gsMatrix& displ); - /// Set the velocity - void setVelocity(gsMatrix& velo); - /// Set the acceleration - void setAcceleration(gsMatrix& accel); - - /// Set time integration method - void setMethod(std::string method); - - /// Toggle verbosity - void verbose() {m_verbose = true; } - - /// Toggle quasi newton method - void quasiNewton() {m_quasiNewton = true; } - - /// Return the displacements - const gsMatrix& displacements() const { return m_displacements;} - /// Return the velocities - const gsMatrix& velocities() const { return m_velocities;} - /// Return the accelerations - const gsMatrix& accelerations() const { return m_accelerations;} - - /// Constructs the solution - void constructSolution(); - - /// Resets the iterations - void resetIteration(); - -protected: - - /// Compute the time-dependent force - gsVector _computeForce(const T & t); - - /// Compute the time-dependent residual - gsVector _computeResidual(const gsVector & U, const T & t); - - /// Compute the time-indepentent jacobian. TODO: Make time-dependent - gsSparseMatrix _computeJacobian(const gsVector & U, const gsVector & dU); - - /// Factorize the matrix \a M - void _factorizeMatrix(const gsSparseMatrix & M); - - /// Solve the system with right-hand side \a F - gsVector _solveSystem(const gsVector & F); - - - /// Initializes some parameters for the class - void initializeCoefficients(); - /// Initializes the solution - void initializeSolution(); - /// Constructs the system for the explicit Euler method - void constructSystemExplEuler(); - /// Constructs the system for the implicit Euler method - void constructSystemImplEuler(); - /// Constructs the block-system for the implicit Euler method - void constructSystemBlock(); - /// Initialize the solution - void initiate(); - - void _step(); - /// Perform a step with the explicit euler method - void stepExplEuler(); - /// Perform a step with the nonlinear explicit euler method - void stepExplEulerNL(); - /// Perform a step with the implicit euler method - void stepImplEuler(); - /// Perform a step with the nonlinear implicit euler method - void stepImplEulerNL(); - /// Perform a step with the nonlinear implicit euler method using block operations - void stepImplEulerNLOp(); - /// Perform a step with the Newmark method - void stepNewmark(); - /// Perform a step with the nonlinear Newmark method - void stepNewmarkNL(); - /// Perform a step with the Bathe method - void stepBathe(); - /// Perform a step with the nonlinear Bathe method - void stepBatheNL(); - /// Perform a step with the Central difference method - void stepCentralDiff(); - /// Perform a step with the nonlinear Central difference method - void stepCentralDiffNL(); - /// Perform a step with the Runge-Kutta 4 method - void stepRK4(); - /// Perform a step with the nonlinear Runge-Kutta 4 method - void stepRK4NL(); - -protected: - - /// Linear solver employed - mutable typename gsSparseSolver::uPtr m_solver; - -protected: - - gsSparseMatrix m_mass; - gsSparseMatrix m_damp; - gsSparseMatrix m_stif; - Jacobian_t m_jacobian; - dJacobian_t m_djacobian; - TResidual_t m_residualFun; - TForce_t m_forceFun; - gsSparseMatrix m_jacMat; - gsMatrix m_resVec; - - size_t m_dofs; - T m_dt; - T m_t; - bool m_NL; - - bool m_verbose; - - bool m_quasiNewton; - - bool m_converged; - - int m_maxIterations; - int m_numIterations; - - bool m_first; - - T m_tolerance; - - gsMatrix m_massInv; - gsMatrix m_forceVec; - - gsSparseMatrix m_sysmat; - gsSparseMatrix m_syseye; - typename gsBlockOp::Ptr m_sysmatBlock; - gsMatrix m_sysvec; - - // For Euler-type methods - gsMatrix m_sol; - gsMatrix m_solOld; - gsMatrix m_dsol; - - // For Newmark-type of methods - gsMatrix m_uNew; - gsMatrix m_vNew; - gsMatrix m_aNew; - - // gsFunction m_dini; - // gsFunction m_vini; - - std::string m_method; - - T m_updateNorm; - T m_residue; - - // void constructSolution(); - gsMatrix m_displacements; - gsMatrix m_velocities; - gsMatrix m_accelerations; - - gsStatus m_status; -}; - - -} // namespace gismo - -////////////////////////////////////////////////// -////////////////////////////////////////////////// - - - -#ifndef GISMO_BUILD_LIB -#include GISMO_HPP_HEADER(gsTimeIntegrator.hpp) -#endif diff --git a/src/gsDynamicSolvers/gsTimeIntegrator.hpp b/src/gsDynamicSolvers/gsTimeIntegrator.hpp deleted file mode 100644 index 88e7b7a..0000000 --- a/src/gsDynamicSolvers/gsTimeIntegrator.hpp +++ /dev/null @@ -1,985 +0,0 @@ -/** @file gsTimeIntegrator.hpp - - @brief Provides temporal solvers for structural analysis problems - - TODO: Split in multiple classes with one base class gsTimeIntegratorBase (as gsALMBase) - TODO: Integrate with XBraid - - This file is part of the G+Smo library. - - This Source Code Form is subject to the terms of the Mozilla Public - License, v. 2.0. If a copy of the MPL was not distributed with this - file, You can obtain one at http://mozilla.org/MPL/2.0/. - - Author(s): H.M. Verhelst (2019-..., TU Delft) -*/ - -#pragma once - -#include -#include - -namespace gismo -{ - template - gsVector gsTimeIntegrator::_computeForce(const T & time) - { - gsVector forceVec; - if (!m_forceFun(time, forceVec)) - throw 2; - return forceVec; - } - - template - gsVector gsTimeIntegrator::_computeResidual(const gsVector & U, const T & time) - { - gsVector resVec; - if (!m_residualFun(U, time, resVec)) - throw 2; - return resVec; - } - - template - gsSparseMatrix gsTimeIntegrator::_computeJacobian(const gsVector & U, const gsVector & deltaU) - { - // Compute Jacobian - gsSparseMatrix m; - if (!m_djacobian(U,deltaU,m)) - throw 2; - return m; - } - - template - void gsTimeIntegrator::_factorizeMatrix(const gsSparseMatrix & M) - { - m_solver->compute(M); - if (m_solver->info()!=gsEigen::ComputationInfo::Success) - { - gsInfo<<"Solver error with code "<info()<<". See Eigen documentation on ComputationInfo \n" - < - gsVector gsTimeIntegrator::_solveSystem(const gsVector & F) - { - try - { - return m_solver->solve(F); - } - catch (...) - { - throw 3; - } - } - - template - void gsTimeIntegrator::initializeCoefficients() - { - m_verbose = false; - m_quasiNewton = false; - m_numIterations = 0; - m_maxIterations = 100; - m_tolerance = 1e-12; - m_converged = false; - m_t = 0; - m_residue = 1.0; - m_updateNorm = 1.0; - m_first = true; - } - - template - void gsTimeIntegrator::initializeSolution() - { - // Defines homogenous solutions as initial solution. Can be overwritten with initialDisplacement() or initialVelocity() later - if (m_method=="ExplEuler" || m_method=="ImplEuler" || m_method=="RK4") - { - m_sol = gsMatrix::Zero(2*m_dofs,1); - m_aNew = gsMatrix::Zero(m_dofs,1); - } - else if ((m_method=="Newmark") || (m_method=="Bathe") || (m_method=="CentralDiff")) - { - m_uNew = gsMatrix::Zero(m_dofs,1); - m_vNew = gsMatrix::Zero(m_dofs,1); - m_aNew = gsMatrix::Zero(m_dofs,1); - } - else - { - gsInfo<<"Method unknown..."; - } - } - - template - void gsTimeIntegrator::resetSolution() - { - if (m_method=="ExplEuler" || m_method=="ImplEuler" || m_method=="RK4") - { - m_sol = gsMatrix::Zero(2*m_dofs,1); - m_aNew = gsMatrix::Zero(m_dofs,1); - } - else if ((m_method=="Newmark") || (m_method=="Bathe") || (m_method=="CentralDiff")) - { - m_uNew = gsMatrix::Zero(m_dofs,1); - m_vNew = gsMatrix::Zero(m_dofs,1); - m_aNew = gsMatrix::Zero(m_dofs,1); - } - - } - - template - void gsTimeIntegrator::setDisplacement(gsMatrix& displ) - { - if (m_method=="ExplEuler" || m_method=="ImplEuler" || m_method=="RK4") - m_sol.block(0,0,m_dofs,1) = displ; - else if ((m_method=="Newmark") || (m_method=="Bathe") || (m_method=="CentralDiff")) - m_uNew = displ; - } - - template - void gsTimeIntegrator::setVelocity(gsMatrix& velo) - { - if (m_method=="ExplEuler" || m_method=="ImplEuler" || m_method=="RK4") - m_sol.block(m_dofs,0,m_dofs,1) = velo; - else if ((m_method=="Newmark") || (m_method=="Bathe") || (m_method=="CentralDiff")) - m_vNew = velo; - } - - template - void gsTimeIntegrator::setAcceleration(gsMatrix& accel) - { - if (m_method=="ExplEuler" || m_method=="ImplEuler" || m_method=="RK4") - // gsWarn<<"Initial acceleration not implemented for method "< - void gsTimeIntegrator::resetIteration() - { - m_numIterations = 0; - m_residue = 1.0; - m_updateNorm = 1.0; - } - - template - void gsTimeIntegrator::constructSystemExplEuler() - { - m_sysmat = gsSparseMatrix(2*m_dofs,2*m_dofs); - m_sysmat.setZero(); - - m_sysvec.setZero(2*m_dofs,1); - - m_massInv = m_mass.toDense().inverse(); - // Explicit Euler Method - gsSparseMatrixImat(m_dofs,m_dofs); - Imat.setIdentity(); - - gsMatrix sub1 = -m_massInv*m_stif; - gsMatrix sub2 = -m_massInv*m_damp; - - // top-right - for (size_t i=0; i < m_dofs; i++) - m_sysmat.insert(i,i+m_dofs) = 1.0; - // bottom-left - for (size_t i=0; i < m_dofs; i++) - { - for (size_t j=0; j < m_dofs; j++) - { - m_sysmat.insert(i+m_dofs,j) = sub1(i,j); - } - } - // bottom-right - for (size_t i=0; i < m_dofs; i++) - { - for (size_t j=0; j < m_dofs; j++) - { - m_sysmat.insert(i+m_dofs,j+m_dofs) = sub2(i,j); - } - } - - m_sysmat *= m_dt; - } - - template - void gsTimeIntegrator::constructSystemImplEuler() - { - constructSystemExplEuler(); - - gsSparseMatrix eye(2*m_dofs,2*m_dofs); - eye.setIdentity(); - m_sysmat = eye-m_sysmat; - - // gsSparseMatrix eye(2*m_dofs,2*m_dofs); - // for (index_t i=0; i < 2*m_dofs; i++) - // m_sysmat(i,i+m_dofs) = 1.0 - m_sysmat(i,i); - - } - - template - void gsTimeIntegrator::constructSystemBlock() - { - m_massInv = m_mass.toDense().inverse(); - - m_sysmatBlock = gsBlockOp<>::make(2,2); - m_sysvec = gsMatrix::Zero(2*m_dofs,1); - - gsSparseMatrix<> eye(m_dofs,m_dofs); - eye.setIdentity(); - // m_sysmatBlock->addOperator(0,0,gsIdentityOp::make(m_dofs) ); - // m_sysmatBlock->addOperator(1,1,gsIdentityOp::make(m_dofs) ); - // m_sysmatBlock->addOperator(0,0,makeMatrixOp(eye) ); - // m_sysmatBlock->addOperator(1,1,makeMatrixOp(eye) ); - m_sysmatBlock->addOperator(0,0,gsIdentityOp::make(m_dofs) ); - m_sysmatBlock->addOperator(0,1,makeMatrixOp(-m_dt*eye) ); - m_sysmatBlock->addOperator(1,0,makeMatrixOp(m_dt*m_massInv*m_stif) ); - m_sysmatBlock->addOperator(1,1,makeMatrixOp(eye + m_dt*m_massInv*m_damp) ); - - gsMatrix uNew = m_sol.topRows(m_dofs); - gsMatrix uOld = m_solOld.topRows(m_dofs); - gsMatrix vNew = m_sol.bottomRows(m_dofs); - gsMatrix vOld = m_solOld.bottomRows(m_dofs); - - m_sysvec.topRows(m_dofs) = uNew - uOld - m_dt*vNew; - m_sysvec.bottomRows(m_dofs) = vNew - vOld + m_dt*(m_massInv*(-m_resVec)); - } - - template - void gsTimeIntegrator::setMethod(std::string string) - { - if ( (string != "ExplEuler") && - (string != "ImplEuler") && - (string != "Newmark") && - (string != "Bathe") && - (string !="CentralDiff") && - (string !="RK4") ) - { - gsInfo<<"Method Unknown... Process terminating\n"; - } - else - { - m_method = string; - initializeSolution(); - } - } - - template - void gsTimeIntegrator::initiate() - { - // BUILD SYSTEM - if (m_method=="ExplEuler") - { - initializeSolution(); - if (!m_NL) - constructSystemExplEuler(); - } - else if (m_method=="ImplEuler") - { - initializeSolution(); - if (!m_NL) - constructSystemImplEuler(); - } - else if (m_method=="Newmark") - { - initializeSolution(); - } - else if (m_method=="Bathe") - { - initializeSolution(); - } - else if (m_method=="CentralDiff") - { - initializeSolution(); - } - else if (m_method=="RK4") - { - initializeSolution(); - } - else - { - gsInfo<<"Method unknown..."; - } - m_solver = gsSparseSolver::get( "LU" ); // todo: make proper options - - // case "ImplEulerBlockOp": - // { - - // } - } - - template - void gsTimeIntegrator::stepExplEuler() - { - // ONLY FOR LINEAR SYSTEM! - m_solOld = m_sol; - gsMatrix uOld = m_solOld.topRows(m_dofs); - gsMatrix vOld = m_solOld.bottomRows(m_dofs); - m_forceVec = _computeForce(m_t); - m_massInv = m_mass.toDense().inverse(); - - m_sol.topRows(m_dofs) += m_dt * vOld; - m_sol.bottomRows(m_dofs) += m_dt * m_massInv * ( m_forceVec - m_stif * uOld - m_damp * vOld); - - m_t += m_dt; - } - - template - void gsTimeIntegrator::stepExplEulerNL() - { - m_solOld = m_sol; - gsMatrix uOld = m_solOld.topRows(m_dofs); - gsMatrix vOld = m_solOld.bottomRows(m_dofs); - m_resVec = _computeResidual(uOld,m_t); - m_massInv = m_mass.toDense().inverse(); - - m_sol.topRows(m_dofs) += m_dt * vOld; - m_sol.bottomRows(m_dofs) += m_dt * m_massInv * ( m_resVec - m_damp * vOld); - - m_t += m_dt; - } - - template - void gsTimeIntegrator::stepImplEuler() - { - this->constructSystemImplEuler(); - m_forceVec = _computeForce(m_t); - m_sysvec.bottomRows(m_dofs) = m_dt*m_massInv*m_forceVec; - - _factorizeMatrix(m_sysmat); - m_sol = _solveSystem(m_sysvec+m_sol); - m_t += m_dt; - } - - template - void gsTimeIntegrator::stepImplEulerNL() - { - gsMatrix solOld = m_sol; - gsMatrix uOld = solOld.topRows(m_dofs); - gsMatrix vOld = solOld.bottomRows(m_dofs); - gsMatrix uNew,vNew; - gsMatrix du = gsMatrix::Zero(m_dofs,1); - - m_sysmat = gsSparseMatrix(2*m_dofs,2*m_dofs); - m_sysmat.setZero(); - - m_sysvec.setZero(2*m_dofs,1); - - gsSparseMatrix eyeN(m_dofs,m_dofs); - eyeN.setIdentity(); - gsSparseMatrix eye2N(2*m_dofs,2*m_dofs); - eye2N.setIdentity(); - - while ( (m_updateNorm>m_tolerance || m_residue>m_tolerance) && m_numIterations<=m_maxIterations ) - { - uNew = m_sol.topRows(m_dofs); - vNew = m_sol.bottomRows(m_dofs); - - if ( (!m_quasiNewton) || (m_numIterations==0) ) - { - du = m_dsol.topRows(m_dofs); - m_jacMat = _computeJacobian(uNew,du); - } - - m_resVec = _computeResidual(uNew,m_t); - - m_sysvec.topRows(m_dofs) = uNew - uOld - m_dt*vNew; - m_sysvec.bottomRows(m_dofs) = m_mass*(vNew - vOld) + m_dt * m_damp * vNew + m_dt*((-m_resVec)); - - gsMatrix sub1 = m_dt*m_jacMat; - gsMatrix sub2 = m_dt*m_damp + m_mass*eyeN; - - m_sysmat.setZero(); - // top-left - for (size_t i=0; i < m_dofs; i++) - m_sysmat.insert(i,i) = 1.0; - // top-right - for (size_t i=0; i < m_dofs; i++) - m_sysmat.insert(i,i+m_dofs) = -m_dt*1.0; - // bottom-left - for (size_t i=0; i < m_dofs; i++) - { - for (size_t j=0; j < m_dofs; j++) - { - m_sysmat.insert(i+m_dofs,j) = sub1(i,j); - } - } - // bottom-right - for (size_t i=0; i < m_dofs; i++) - { - for (size_t j=0; j < m_dofs; j++) - { - m_sysmat.insert(i+m_dofs,j+m_dofs) = sub2(i,j); - } - } - - m_residue = m_sysvec.norm(); - _factorizeMatrix(m_sysmat); - m_dsol = _solveSystem(-m_sysvec); - - m_sol += m_dsol; - m_updateNorm = m_dsol.norm() / m_sol.norm(); - - if (m_verbose) - { - gsInfo <<"\tIteration: " << m_numIterations - <<", residue: " << m_residue - <<", update norm: " << m_updateNorm - <<"\n"; - } - m_numIterations++; - } - - m_converged = true; - if ((m_numIterations==m_maxIterations+1) && m_verbose) - { - gsInfo<<"\tWARNING: Maximum number of iterations reached"<<"\n"; - m_converged = false; - } - - m_t += m_dt; - - resetIteration(); - } - - template - void gsTimeIntegrator::stepImplEulerNLOp() - { - gsMatrix solOld = m_sol; - gsMatrix dsol; - gsMatrix uOld = solOld.topRows(m_dofs); - gsMatrix vOld = solOld.bottomRows(m_dofs); - gsMatrix uNew,vNew; - gsMatrix du; - m_dsol = gsMatrix::Zero(2*m_dofs,1); - - gsBlockOp<>::Ptr Amat=gsBlockOp<>::make(2,2); - - gsSparseMatrix eye(m_dofs,m_dofs); - eye.setIdentity(); - - // m_sysmatBlock=gsBlockOp<>::make(2,2); - m_sysvec = gsMatrix::Zero(2*m_dofs,1); - - gsGMRes<> gmres(Amat); - gmres.setMaxIterations(100); - gmres.setTolerance(1e-8); - while ( (m_updateNorm>m_tolerance || m_residue>m_tolerance) && m_numIterations<=m_maxIterations ) - { - uNew = m_sol.topRows(m_dofs); - vNew = m_sol.bottomRows(m_dofs); - - if ( (!m_quasiNewton) || (m_numIterations==0) ) - { - du = m_dsol.topRows(m_dofs); - m_jacMat = _computeJacobian(uNew,du); - } - - m_resVec = _computeResidual(uNew,m_t); - - m_sysvec.topRows(m_dofs) = uNew - uOld - m_dt*vNew; - m_sysvec.bottomRows(m_dofs) = m_mass*(vNew - vOld) + m_dt * m_damp * vNew + m_dt*(-m_resVec); - - // m_jacMat = _computeJacobian(uNew,du); - - Amat->addOperator(0,0,gsIdentityOp::make(m_dofs) ); - Amat->addOperator(0,1,makeMatrixOp(-m_dt*eye) ); - Amat->addOperator(1,0,makeMatrixOp(m_dt*m_jacMat) ); - Amat->addOperator(1,1,makeMatrixOp(m_mass + m_dt*m_damp) ); - - m_residue= m_sysvec.norm(); - gmres.solve(-m_sysvec,dsol); - - m_dsol = dsol; - - m_sol += m_dsol; - m_updateNorm = dsol.norm() / m_sol.norm(); - - if (m_verbose) - { - gsInfo <<"\tIteration: " << m_numIterations - <<", residue: " << m_residue - <<", update norm: " << m_updateNorm - <<"\n"; - } - m_numIterations++; - } - - m_converged = true; - if ((m_numIterations==m_maxIterations+1) && m_verbose) - { - gsInfo<<"\tWARNING: Maximum number of iterations reached"<<"\n"; - m_converged = false; - } - - m_t += m_dt; - - resetIteration(); - } - - template - void gsTimeIntegrator::stepNewmark() - { - T gamma = 1; - m_t += gamma*m_dt; - - gsMatrix uOld,vOld,aOld; - uOld = m_uNew; - vOld = m_vNew; - aOld = m_aNew; - - m_forceVec = _computeForce(m_t); - - gsSparseMatrix<> lhs = m_stif + 4/(math::pow(gamma*m_dt,2))*m_mass + 2/(gamma*m_dt)*m_damp; - gsMatrix<> rhs = m_forceVec + m_mass*(4/(math::pow(gamma*m_dt,2))*(uOld)+4/(gamma*m_dt)*vOld+aOld) + m_damp*(2/(gamma*m_dt)*(uOld)+vOld); - - _factorizeMatrix(lhs); - m_sol = _solveSystem(rhs); - - m_uNew = m_sol; - - m_vNew = 2/(gamma*m_dt)*(m_uNew-uOld) - vOld; - m_aNew = (m_uNew-uOld-vOld*m_dt)*4/(math::pow(gamma*m_dt,2)) - aOld; - } - - template - void gsTimeIntegrator::stepNewmarkNL() - { - T gamma = 1; - m_t += gamma*m_dt; - - gsMatrix uOld,vOld,aOld; - uOld = m_uNew; - vOld = m_vNew; - aOld = m_aNew; - gsMatrix du; - m_dsol = gsMatrix::Zero(m_dofs,1); - - - gsMatrix<> rhs; - gsSparseMatrix<> lhs; - - m_resVec = _computeResidual(m_uNew,m_t); - rhs = m_resVec - m_mass*(4/(math::pow(gamma*m_dt,2))*(m_uNew-uOld)-4/(gamma*m_dt)*vOld-aOld) - m_damp*(2/(gamma*m_dt)*(m_uNew-uOld)-vOld); - T res0 = rhs.norm(); - - while ( (m_updateNorm>m_tolerance || m_residue/res0>m_tolerance) && m_numIterations<=m_maxIterations ) - { - if ( (!m_quasiNewton) || (m_numIterations==0) ) - { - du = m_dsol.topRows(m_dofs); - m_jacMat = _computeJacobian(m_uNew,du); - } - - lhs = m_jacMat + 4/(math::pow(gamma*m_dt,2))*m_mass + 2/(gamma*m_dt)*m_damp; - _factorizeMatrix(lhs); - m_dsol = _solveSystem(rhs); - m_updateNorm = m_dsol.norm()/m_uNew.norm(); - m_uNew += m_dsol; - - m_resVec = _computeResidual(m_uNew,m_t); - rhs = m_resVec - m_mass*(4/(math::pow(gamma*m_dt,2))*(m_uNew-uOld)-4/(gamma*m_dt)*vOld-aOld) - m_damp*(2/(gamma*m_dt)*(m_uNew-uOld)-vOld); - m_residue = rhs.norm(); - - if (m_verbose) - { - gsInfo <<"\tIteration: " << m_numIterations - <<", residue: " << m_residue - <<", update norm: " << m_updateNorm - <<"\n"; - } - m_numIterations++; - } - // Construct velocities and accelerations - m_vNew = 2/(gamma*m_dt)*(m_uNew-uOld) - vOld; - m_aNew = (m_uNew-uOld-vOld*m_dt)*4/(math::pow(gamma*m_dt,2)) - aOld; - - m_converged = true; - if ((m_numIterations==m_maxIterations+1) && m_verbose) - { - gsInfo<<"\tWARNING: Maximum number of iterations reached"<<"\n"; - m_converged = false; - } - - resetIteration(); - } - - template - void gsTimeIntegrator::stepBathe() - { - T gamma = 0.5; - m_t += gamma*m_dt; - - gsMatrix uOld,vOld,aOld; - gsMatrix uStep,vStep,aStep; - uOld = m_uNew; - vOld = m_vNew; - aOld = m_aNew; - - m_forceVec = _computeForce(m_t); - - gsSparseMatrix<> lhs = m_stif + 4/(math::pow(gamma*m_dt,2))*m_mass + 2/(gamma*m_dt)*m_damp; - gsMatrix<> rhs = m_forceVec + m_mass*(4/(math::pow(gamma*m_dt,2))*(uOld)+4/(gamma*m_dt)*vOld+aOld) + m_damp*(2/(gamma*m_dt)*(uOld)+vOld);; - _factorizeMatrix(lhs); - m_sol = _solveSystem(rhs); - m_uNew = m_sol; - - m_vNew = 2/(gamma*m_dt)*(m_uNew-uOld) - vOld; - m_aNew = (m_uNew-uOld-vOld*m_dt)*4/(math::pow(gamma*m_dt,2)) - aOld; - - uStep = m_uNew; - vStep = m_vNew; - aStep = m_aNew; - - T c1 = (1-gamma)/(gamma*m_dt); - T c2 = (-1)/((1-gamma)*gamma*m_dt); - T c3 = (2-gamma)/((1-gamma)*m_dt); - - m_t += (1-gamma)*m_dt; - - lhs = m_stif + c3*c3*m_mass + c3*m_damp; - rhs = m_forceVec - m_mass*(c1*vOld+c2*vStep+c1*c3*uOld+c3*c2*uStep) - m_damp*(c2*uStep+c1*uOld); - _factorizeMatrix(lhs); - m_sol = _solveSystem(rhs); - m_uNew = m_sol; - m_residue = rhs.norm(); - - m_vNew = c1*uOld + c2*uStep + c3*m_uNew; - m_aNew = c1*vOld + c2*vStep + c3*m_vNew; - - } - - template - void gsTimeIntegrator::stepBatheNL() - { - T gamma = 0.5; - m_t += gamma*m_dt; - - gsMatrix uOld,vOld,aOld; - gsMatrix uStep,vStep,aStep; - uOld = m_uNew; - vOld = m_vNew; - aOld = m_aNew; - gsMatrix du; - m_dsol = gsMatrix::Zero(m_dofs,1); - - if (m_verbose) - { - gsInfo << "\tStage 1:\n"; - } - while ( (m_updateNorm>m_tolerance || m_residue>m_tolerance) && m_numIterations<=m_maxIterations ) - { - if ( (!m_quasiNewton) || (m_numIterations==0) ) - { - du = m_dsol.topRows(m_dofs); - m_jacMat = _computeJacobian(m_uNew,du); - } - - m_resVec = _computeResidual(m_uNew,m_t); - - gsSparseMatrix<> lhs = m_jacMat + 4/(math::pow(gamma*m_dt,2))*m_mass + 2/(gamma*m_dt)*m_damp; - gsMatrix<> rhs = m_resVec - m_mass*(4/(math::pow(gamma*m_dt,2))*(m_uNew-uOld)-4/(gamma*m_dt)*vOld-aOld) - m_damp*(2/(gamma*m_dt)*(m_uNew-uOld)-vOld); - _factorizeMatrix(lhs); - m_dsol = _solveSystem(rhs); - m_uNew += m_dsol; - m_updateNorm = m_dsol.norm() / m_uNew.norm(); - m_residue = rhs.norm(); - - if (m_verbose) - { - gsInfo <<"\tIteration: " << m_numIterations - <<", residue: " << m_residue - <<", update norm: " << m_updateNorm - <<"\n"; - } - m_numIterations++; - } - // Construct velocities and accelerations - m_vNew = 2/(gamma*m_dt)*(m_uNew-uOld) - vOld; - m_aNew = (m_uNew-uOld-vOld*m_dt)*4/(math::pow(gamma*m_dt,2)) - aOld; - uStep = m_uNew; - vStep = m_vNew; - aStep = m_aNew; - - T c1 = (1-gamma)/(gamma*m_dt); - T c2 = (-1)/((1-gamma)*gamma*m_dt); - T c3 = (2-gamma)/((1-gamma)*m_dt); - - if (m_verbose) - { - gsInfo << "\tStage 2:\n"; - } - - m_t += (1-gamma)*m_dt; - resetIteration(); - - while ( (m_updateNorm>m_tolerance || m_residue>m_tolerance) && m_numIterations<=m_maxIterations ) - { - if ( (!m_quasiNewton) || (m_numIterations==0) ) - { - du = m_dsol.topRows(m_dofs); - m_jacMat = _computeJacobian(m_uNew,du); - } - - m_resVec = _computeResidual(m_uNew,m_t); - - gsSparseMatrix lhs = m_jacMat + c3*c3*m_mass + c3*m_damp; - gsMatrix rhs = m_resVec - m_mass*(c1*vOld+c2*vStep+c1*c3*uOld+c3*c2*uStep+c3*c3*m_uNew) - m_damp*(c1*uOld+c2*uStep+c3*m_uNew); - _factorizeMatrix(lhs); - m_dsol = _solveSystem(rhs); - m_uNew += m_dsol; - m_updateNorm = m_dsol.norm() / m_uNew.norm(); - m_residue = rhs.norm(); - - if (m_verbose) - { - gsInfo <<"\tIteration: " << m_numIterations - <<", residue: " << m_residue - <<", update norm: " << m_updateNorm - <<"\n"; - } - m_numIterations++; - } - - m_vNew = c1*uOld + c2*uStep + c3*m_uNew; - m_aNew = c1*vOld + c2*vStep + c3*m_vNew; - - m_converged = true; - if ((m_numIterations==m_maxIterations+1) && m_verbose) - { - gsInfo<<"\tWARNING: Maximum number of iterations reached"<<"\n"; - m_converged = false; - } - - resetIteration(); - } - - template - void gsTimeIntegrator::stepCentralDiff() - { - gsMatrix uOld,vOld,aOld; - m_massInv = m_mass.toDense().inverse(); - - uOld = m_uNew; // u_n - vOld = m_vNew; // v_n+1/2 - aOld = m_aNew; // a_n - m_forceVec = _computeForce(m_t); - - if (m_first) - { - aOld = m_massInv * (m_forceVec - m_stif * uOld - m_damp * vOld); - vOld = vOld + m_dt / 2 * aOld; // v_1/2 = v_0 + dt/2*a_0 - m_first = false; - } - - m_uNew = uOld + m_dt * vOld; - m_aNew = m_massInv * (m_forceVec - m_stif * m_uNew - m_damp * vOld); - m_vNew = vOld + m_dt * m_aNew; // v_n+1/2 = v_n-1/2 + dt*a_n - - m_t+=m_dt; - } - - template - void gsTimeIntegrator::stepCentralDiffNL() - { - gsMatrix uOld,vOld,aOld; - m_massInv = m_mass.toDense().inverse(); - - uOld = m_uNew; // u_n - vOld = m_vNew; // v_n+1/2 - aOld = m_aNew; // a_n - - if (m_first) - { - m_resVec = _computeResidual(uOld,m_t); - aOld = m_massInv * (m_resVec - m_damp * vOld); - vOld = vOld + m_dt / 2 * aOld; // v_1/2 = v_0 + dt/2*a_0 - m_first = false; - } - - m_uNew = uOld + m_dt * vOld; - m_t+=m_dt; - m_resVec = _computeResidual(uOld,m_t); - m_aNew = m_massInv * (m_resVec - m_damp * vOld); - m_vNew = vOld + m_dt * m_aNew; // v_n+1/2 = v_n-1/2 + dt*a_n - - } - - template - void gsTimeIntegrator::stepRK4() - { - m_massInv = m_mass.toDense().inverse(); - gsVector k1(2*m_dofs), k2(2*m_dofs), k3(2*m_dofs), k4(2*m_dofs); - gsVector uTmp, vTmp; - m_solOld = m_sol; - - gsVector uOld = m_solOld.topRows(m_dofs); - gsVector vOld = m_solOld.bottomRows(m_dofs); - - m_resVec = _computeForce(m_t) - m_stif * uOld; - k1.topRows(m_dofs) = vOld; - k1.bottomRows(m_dofs) = m_massInv * ( m_resVec - m_damp * vOld); - - uTmp = uOld + m_dt/2. * k1.topRows(m_dofs); - vTmp = vOld + m_dt/2. * k1.bottomRows(m_dofs); - m_resVec = _computeForce(m_t + m_dt/2.) - m_stif * uTmp; - k2.topRows(m_dofs) = vTmp; - k2.bottomRows(m_dofs) = m_massInv * ( m_resVec - m_damp * vTmp); - - uTmp = uOld + m_dt/2. * k2.topRows(m_dofs); - vTmp = vOld + m_dt/2. * k2.bottomRows(m_dofs); - m_resVec = _computeForce(m_t + m_dt/2.) - m_stif * uTmp; - k3.topRows(m_dofs) = vTmp; - k3.bottomRows(m_dofs) = m_massInv * ( m_resVec - m_damp * vTmp); - - uTmp = uOld + m_dt * k3.topRows(m_dofs); - vTmp = vOld + m_dt * k3.bottomRows(m_dofs); - m_resVec = _computeForce(m_t + m_dt) - m_stif * uTmp; - k4.topRows(m_dofs) = vTmp; - k4.bottomRows(m_dofs) = m_massInv * ( m_resVec - m_damp * vTmp); - - m_sol += 1./6. * m_dt * ( k1 + 2.*k2 + 2.*k3 + k4 ); - m_t += m_dt; - } - - template - void gsTimeIntegrator::stepRK4NL() - { - m_massInv = m_mass.toDense().inverse(); - gsVector k1(2*m_dofs), k2(2*m_dofs), k3(2*m_dofs), k4(2*m_dofs); - gsVector uTmp, vTmp; - m_solOld = m_sol; - - gsVector uOld = m_solOld.topRows(m_dofs); - gsVector vOld = m_solOld.bottomRows(m_dofs); - - m_resVec = _computeResidual(uOld,m_t); - k1.topRows(m_dofs) = vOld; - k1.bottomRows(m_dofs) = m_massInv * ( m_resVec - m_damp * vOld); - - uTmp = uOld + m_dt/2. * k1.topRows(m_dofs); - vTmp = vOld + m_dt/2. * k1.bottomRows(m_dofs); - m_resVec = _computeResidual(uTmp,m_t + m_dt/2.); - k2.topRows(m_dofs) = vTmp; - k2.bottomRows(m_dofs) = m_massInv * ( m_resVec - m_damp * vTmp); - - uTmp = uOld + m_dt/2. * k2.topRows(m_dofs); - vTmp = vOld + m_dt/2. * k2.bottomRows(m_dofs); - m_resVec = _computeResidual(uTmp,m_t + m_dt/2.); - k3.topRows(m_dofs) = vTmp; - k3.bottomRows(m_dofs) = m_massInv * ( m_resVec - m_damp * vTmp); - - uTmp = uOld + m_dt * k3.topRows(m_dofs); - vTmp = vOld + m_dt * k3.bottomRows(m_dofs); - m_resVec = _computeResidual(uTmp,m_t + m_dt); - k4.topRows(m_dofs) = vTmp; - k4.bottomRows(m_dofs) = m_massInv * ( m_resVec - m_damp * vTmp); - - m_sol += 1./6. * m_dt * ( k1 + 2.*k2 + 2.*k3 + k4 ); - - m_t += m_dt; - } - - template - gsStatus gsTimeIntegrator::step() - { - try - { - _step(); - m_status = gsStatus::Success; - } - catch (int errorCode) - { - if (errorCode==1) - m_status = gsStatus::NotConverged; - else if (errorCode==2) - m_status = gsStatus::AssemblyError; - else if (errorCode==3) - m_status = gsStatus::SolverError; - else - m_status = gsStatus::OtherError; - } - catch (...) - { - m_status = gsStatus::OtherError; - } - return m_status; - } - - template - void gsTimeIntegrator::_step() - { - if (m_verbose) - gsInfo << "time = "< - void gsTimeIntegrator::constructSolution() - { - if ((m_method=="ExplEuler") || (m_method=="ImplEuler") || (m_method=="RK4")) - { - m_displacements = m_sol.topRows(m_dofs); - m_velocities = m_sol.bottomRows(m_dofs); - m_accelerations = m_aNew; - } - else if ((m_method=="Newmark") || (m_method=="Bathe") || (m_method=="CentralDiff")) - { - m_displacements = m_uNew; - m_velocities = m_vNew; - m_accelerations = m_aNew; - } - } - -} // namespace gismo \ No newline at end of file diff --git a/src/gsDynamicSolvers/gsTimeIntegrator_.cpp b/src/gsDynamicSolvers/gsTimeIntegrator_.cpp deleted file mode 100644 index f284352..0000000 --- a/src/gsDynamicSolvers/gsTimeIntegrator_.cpp +++ /dev/null @@ -1,9 +0,0 @@ -#include - -#include -#include - -namespace gismo -{ - CLASS_TEMPLATE_INST gsTimeIntegrator; -} diff --git a/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h b/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h index 77132a2..f571432 100644 --- a/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h +++ b/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h @@ -16,47 +16,75 @@ namespace gismo { - enum struct gsStatus - { - Success, ///< Successful - NotConverged, ///< Step did not converge - AssemblyError, ///< Assembly problem in step - SolverError, ///< Assembly problem in step - NotStarted, ///< ALM has not started yet - OtherError ///< Other error - }; - - template - struct gsStructuralAnalysisOps - { - /// Force - typedef std::function < bool ( gsVector & )> Force_t; - /// Time-dependent force - typedef std::function < bool ( const T, gsVector & )> TForce_t; - - /// Residual, Fint-Fext - typedef std::function < bool ( gsVector const &, gsVector & )> Residual_t; - /// Arc-Length Residual, Fint-lambda*Fext - typedef std::function < bool ( gsVector const &, const T, gsVector & )> ALResidual_t; - /// Time-dependent Residual Fint(t)-Fext(t) - typedef std::function < bool ( gsVector const &, const T, gsVector & )> TResidual_t; - - /// Mass matrix - typedef std::function < bool ( gsSparseMatrix & ) > Mass_t; - /// Time-dependent mass matrix - typedef std::function < bool ( const T, gsSparseMatrix & ) > TMass_t; - /// Damping matrix - typedef std::function < bool ( gsVector const &, gsSparseMatrix & ) > Damping_t; - /// Time-dependent Damping matrix - typedef std::function < bool ( gsVector const &, const T, gsSparseMatrix & ) > TDamping_t; - - /// Stiffness matrix - typedef std::function < bool ( gsSparseMatrix & ) > Stiffness_t; - /// Jacobian - typedef std::function < bool ( gsVector const &, gsSparseMatrix & ) > Jacobian_t; - /// Jacobian - typedef std::function < bool ( gsVector const &, const T, gsSparseMatrix & ) > TJacobian_t; - /// Jacobian with solution update as argument - typedef std::function < bool ( gsVector const &, gsVector const &, gsSparseMatrix & ) > dJacobian_t; - }; -} // namespace gismo + +enum struct gsStatus +{ + Success, ///< Successful + NotConverged, ///< Step did not converge + AssemblyError, ///< Assembly problem in step + SolverError, ///< Assembly problem in step + NotStarted, ///< ALM has not started yet + OtherError ///< Other error +}; + +// ALTERNATIVE IMPLEMENTATION USING FUNCTORS +// template +// struct DynamicFunctor +// { +// Force(); + +// gsVector& operator(T time = 0) {}; + +// virtual index_t rows() = 0; +// virtual index_t cols() = 0; +// }; + + +// template +// struct DynamicForce : public DynamicFunctor +// { +// DynamicForce(... ) {} + + +// }; + +/** + * @brief Operators for the gsStructuralAnalysis module + * + * @tparam T Double type + */ +template +struct gsStructuralAnalysisOps +{ + /// Force + typedef std::function < bool ( gsVector & )> Force_t; + /// Time-dependent force + typedef std::function < bool ( const T, gsVector & )> TForce_t; + + /// Residual, Fint-Fext + typedef std::function < bool ( gsVector const &, gsVector & )> Residual_t; + /// Arc-Length Residual, Fint-lambda*Fext + typedef std::function < bool ( gsVector const &, const T, gsVector & )> ALResidual_t; + /// Time-dependent Residual Fint(t)-Fext(t) + typedef std::function < bool ( gsVector const &, const T, gsVector & )> TResidual_t; + + /// Mass matrix + typedef std::function < bool ( gsSparseMatrix & ) > Mass_t; + /// Time-dependent mass matrix + typedef std::function < bool ( const T, gsSparseMatrix & ) > TMass_t; + /// Damping matrix + typedef std::function < bool ( gsVector const &, gsSparseMatrix & ) > Damping_t; + /// Time-dependent Damping matrix + typedef std::function < bool ( gsVector const &, const T, gsSparseMatrix & ) > TDamping_t; + + /// Stiffness matrix + typedef std::function < bool ( gsSparseMatrix & ) > Stiffness_t; + /// Jacobian + typedef std::function < bool ( gsVector const &, gsSparseMatrix & ) > Jacobian_t; + /// Jacobian + typedef std::function < bool ( gsVector const &, const T, gsSparseMatrix & ) > TJacobian_t; + /// Jacobian with solution update as argument + typedef std::function < bool ( gsVector const &, gsVector const &, gsSparseMatrix & ) > dJacobian_t; +}; + +} diff --git a/tutorials/CMakeLists.txt b/tutorials/CMakeLists.txt new file mode 100644 index 0000000..54267b8 --- /dev/null +++ b/tutorials/CMakeLists.txt @@ -0,0 +1,26 @@ +###################################################################### +## CMakeLists.txt --- gsStructuralAnalysis/tutorials +## This file is part of the G+Smo library. +## +## Author: Angelos Mantzaflaris, Hugo Verhelst +## Copyright (C) 2023 +###################################################################### + +# add tutorial files +aux_cpp_directory(${CMAKE_CURRENT_SOURCE_DIR} FILES) +foreach(file ${FILES}) + add_gismo_executable(${file}) + get_filename_component(tarname ${file} NAME_WE) # name without extension + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}") + if(GISMO_BUILD_EXAMPLES) + set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") + else(GISMO_BUILD_EXAMPLES) + set_target_properties(${tarname} PROPERTIES + FOLDER "${PROJECT_NAME}" + EXCLUDE_FROM_ALL TRUE) + endif(GISMO_BUILD_EXAMPLES) + add_dependencies(${PROJECT_NAME}-tutorials ${tarname}) + add_dependencies(${PROJECT_NAME}-all ${tarname}) + # install the example executables (optionally) + install(TARGETS ${tarname} DESTINATION "${BIN_INSTALL_DIR}" COMPONENT exe OPTIONAL) +endforeach(file ${FILES}) \ No newline at end of file diff --git a/tutorials/nonlinear_shell_arcLength.cpp b/tutorials/nonlinear_shell_arcLength.cpp new file mode 100644 index 0000000..f05de5b --- /dev/null +++ b/tutorials/nonlinear_shell_arcLength.cpp @@ -0,0 +1,189 @@ +/** @file gsStructuralAnalysis/tutorials/nonlinear_shell_arclength.cpp + + @brief Tutorial for assembling the Kirchhoff-Love shell + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M.Verhelst +*/ + +#include + +#ifdef gsKLShell_ENABLED +#include +#endif + +#include +#include + + +using namespace gismo; + +// Choose among various shell examples, default = Thin Plate +int main(int argc, char *argv[]) +{ +#ifdef gsKLShell_ENABLED + //! [Parse command line] + bool plot = false; + index_t numRefine = 1; + index_t numElevate = 1; + + gsCmdLine cmd("Linear shell tutorial."); + cmd.addInt( "e", "degreeElevation","Number of degree elevation steps to perform before solving (0: equalize degree in all directions)", numElevate ); + cmd.addInt( "r", "uniformRefine", "Number of Uniform h-refinement steps to perform before solving", numRefine ); + cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + //! [Parse command line] + + //! [Read geometry] + // Initialize [ori]ginal and [def]ormed geometry + gsMultiPatch<> ori; + gsReadFile<>("surfaces/paraboloid.xml", ori); + ori.computeTopology(); + //! [Read geometry] + + //! [Initialize geometry] + // p-refine + if (numElevate!=0) + ori.degreeElevate(numElevate); + // h-refine + for (int r =0; r < numRefine; ++r) + ori.uniformRefine(); + //! [Initialize geometry] + + //! [Construct basis] + gsMultiBasis<> bases(ori); + //! [Construct basis] + + //! [Set boundary conditions] + // Define the boundary conditions object + gsBoundaryConditions<> bc; + // Set the geometry map for computation of the Dirichlet BCs + bc.setGeoMap(ori); + + // Set the boundary conditions + bc.addCornerValue(boundary::southwest, 0.0, 0, -1); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::southeast, 0.0, 0, -1); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northwest, 0.0, 0, -1); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northeast, 0.0, 0, -1); // (corner,value, patch, unknown) + + gsPointLoads pLoads = gsPointLoads(); + gsVector<> point(2); + gsVector<> load (3); + point<< 0.5, 0.5 ; load << 0, 0.0, -1e4 ; + pLoads.addLoad(point, load, 0 ); + //! [Set boundary conditions] + + //! [Set surface force] + // The surface force is defined in the physical space, i.e. 3D + gsFunctionExpr<> force("0","0","0",3); + //! [Set surface force] + + + //! [Define the material matrix class] + // Define material parameters + // The material parameters are defined in the physical domain as well (!) + gsConstantFunction<> E (1.E9,3); + gsConstantFunction<> nu(0.45,3); + gsConstantFunction<> t (1E-2,3); + + // Define a linear material, see \ref gsMaterialMatrixLinear.h + // The first parameter is the physical domain dimension + + // gsMaterialMatrixLinear<3,real_t> materialMatrix(ori,t,E,nu); + + std::vector*> parameters{&E,&nu}; + gsOptionList options; + options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",1); + options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",1); + gsMaterialMatrixBase * materialMatrix = getMaterialMatrix<3,real_t>(ori,t,parameters,options); + //! [Define the material matrix class] + + //! [Define the assembler] + // Define the assembler. + // The first parameter is the physical domain dimension + // The second parameter is the real type + // The third parameter controls the bending stiffness + gsThinShellAssembler<3, real_t, true > assembler(ori,bases,bc,force,materialMatrix); + assembler.setPointLoads(pLoads); + //! [Define the assembler] + + //! [Assemble linear part] + assembler.assemble(); + gsSparseMatrix<> K = assembler.matrix(); + gsVector<> F = assembler.rhs(); + //! [Assemble linear part] + + //! [Define nonlinear jacobian] + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler](gsVector const &x, gsSparseMatrix & m) + { + gsMultiPatch<> def; + assembler.constructSolution(x,def); + ThinShellAssemblerStatus status = assembler.assembleMatrix(def); + m = assembler.matrix(); + return status==ThinShellAssemblerStatus::Success; + }; + //! [Define nonlinear jacobian] + + //! [Define nonlinear residual functions] + // Function for the Residual + gsStructuralAnalysisOps::ALResidual_t ALResidual = [&assembler,&F](gsVector const &x, real_t lam, gsVector & v) + { + gsMultiPatch<> def; + assembler.constructSolution(x,def); + ThinShellAssemblerStatus status = assembler.assembleVector(def); + v = F - lam * F - assembler.rhs(); // assembler rhs - force = Finternal + return status == ThinShellAssemblerStatus::Success; + }; + //! [Define nonlinear residual functions] + + //! [Set ALM solver] + gsInfo<<"Solving system with "< solver(Jacobian,ALResidual,F); + solver.options().setSwitch("Verbose",true); + solver.setLength(0.01); + solver.applyOptions(); + solver.initialize(); + //! [Set ALM solver] + + //! [Solve nonlinear problem] + index_t step = 50; + gsParaviewCollection collection("Deformation"); + gsVector<> solVector; + gsMultiPatch<> displ, def; + for (index_t k=0; k solField(def, displ); + std::string outputName = "Deformation" + std::to_string(k) + "_"; + gsWriteParaview<>( solField, outputName, 1000, true); + collection.addPart(outputName + "0.vts",k); + } + } + + //! [Save the paraview collection] + if (plot) + collection.save(); + //! [Save the paraview collection] + + return EXIT_SUCCESS; +#else + GISMO_ERROR("The tutorial needs to be compiled with gsKLShell enabled"); + return EXIT_FAILED; +#endif +}// end main \ No newline at end of file diff --git a/tutorials/nonlinear_shell_dynamic.cpp b/tutorials/nonlinear_shell_dynamic.cpp new file mode 100644 index 0000000..aa317c5 --- /dev/null +++ b/tutorials/nonlinear_shell_dynamic.cpp @@ -0,0 +1,211 @@ +/** @file gsKLShell/tutorials/linear_shell.cpp + + @brief Tutorial for assembling the Kirchhoff-Love shell + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M.Verhelst +*/ + +#include + +#ifdef gsKLShell_ENABLED +#include +#endif + +#include +#include + + +using namespace gismo; + +// Choose among various shell examples, default = Thin Plate +int main(int argc, char *argv[]) +{ +#ifdef gsKLShell_ENABLED + //! [Parse command line] + bool plot = false; + index_t numRefine = 1; + index_t numElevate = 1; + + gsCmdLine cmd("Linear shell tutorial."); + cmd.addInt( "e", "degreeElevation","Number of degree elevation steps to perform before solving (0: equalize degree in all directions)", numElevate ); + cmd.addInt( "r", "uniformRefine", "Number of Uniform h-refinement steps to perform before solving", numRefine ); + cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + //! [Parse command line] + + //! [Read geometry] + // Initialize [ori]ginal and [def]ormed geometry + gsMultiPatch<> ori; + gsReadFile<>("surfaces/paraboloid.xml", ori); + ori.computeTopology(); + //! [Read geometry] + + //! [Initialize geometry] + // p-refine + if (numElevate!=0) + ori.degreeElevate(numElevate); + // h-refine + for (int r =0; r < numRefine; ++r) + ori.uniformRefine(); + //! [Initialize geometry] + + //! [Construct basis] + gsMultiBasis<> bases(ori); + //! [Construct basis] + + //! [Set boundary conditions] + // Define the boundary conditions object + gsBoundaryConditions<> bc; + // Set the geometry map for computation of the Dirichlet BCs + bc.setGeoMap(ori); + + // Set the boundary conditions + bc.addCornerValue(boundary::southwest, 0.0, 0, -1); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::southeast, 0.0, 0, -1); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northwest, 0.0, 0, -1); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northeast, 0.0, 0, -1); // (corner,value, patch, unknown) + + gsPointLoads pLoads = gsPointLoads(); + gsVector<> point(2); + gsVector<> load (3); + point<< 0.5, 0.5 ; load << 0, 0.0, -1e4 ; + pLoads.addLoad(point, load, 0 ); + //! [Set boundary conditions] + + //! [Set surface force] + // The surface force is defined in the physical space, i.e. 3D + gsFunctionExpr<> force("0","0","0",3); + //! [Set surface force] + + + //! [Define the material matrix class] + // Define material parameters + // The material parameters are defined in the physical domain as well (!) + gsConstantFunction<> E (1.E9,3); + gsConstantFunction<> nu (0.45,3); + gsConstantFunction<> t (1E-2,3); + gsConstantFunction<> rho(1.E5,3); + + // Define a linear material, see \ref gsMaterialMatrixLinear.h + // The first parameter is the physical domain dimension + + // gsMaterialMatrixLinear<3,real_t> materialMatrix(ori,t,E,nu); + + std::vector*> parameters{&E,&nu}; + gsOptionList options; + options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",1); + options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",1); + gsMaterialMatrixBase * materialMatrix = getMaterialMatrix<3,real_t>(ori,t,parameters,rho,options); + //! [Define the material matrix class] + + //! [Define the assembler] + // Define the assembler. + // The first parameter is the physical domain dimension + // The second parameter is the real type + // The third parameter controls the bending stiffness + gsThinShellAssembler<3, real_t, true > assembler(ori,bases,bc,force,materialMatrix); + assembler.setPointLoads(pLoads); + //! [Define the assembler] + + //! [Assemble linear part] + assembler.assemble(); + gsSparseMatrix<> K = assembler.matrix(); + gsVector<> F = assembler.rhs(); + + assembler.assembleMass(); + gsSparseMatrix<> M = assembler.matrix(); + //! [Assemble linear part] + + //! [Define nonlinear residual functions] + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler](gsVector const &x, gsSparseMatrix & m) + { + gsMultiPatch<> def; + assembler.constructSolution(x,def); + ThinShellAssemblerStatus status = assembler.assembleMatrix(def); + m = assembler.matrix(); + return status==ThinShellAssemblerStatus::Success; + }; + // Function for the Residual + gsStructuralAnalysisOps::Residual_t Residual = [&assembler](gsVector const &x, gsVector & v) + { + gsMultiPatch<> def; + assembler.constructSolution(x,def); + ThinShellAssemblerStatus status = assembler.assembleVector(def); + v = assembler.rhs(); + return status==ThinShellAssemblerStatus::Success; + }; + //! [Define nonlinear residual functions] + + //! [Define damping and mass matrices] + gsSparseMatrix<> C = gsSparseMatrix<>(assembler.numDofs(),assembler.numDofs()); + gsStructuralAnalysisOps::Damping_t Damping = [&C](const gsVector &, gsSparseMatrix & m) { m = C; return true; }; + gsStructuralAnalysisOps::Mass_t Mass = [&M]( gsSparseMatrix & m) { m = M; return true; }; + //! [Define damping and mass matrices] + + + //! [Set dynamic solver] + gsInfo<<"Solving system with "< solver(Mass,Damping,Jacobian,Residual); + solver.options().setSwitch("Verbose",true); + //! [Set dynamic solver] + + //! [Initialize solution] + size_t N = assembler.numDofs(); + gsVector<> U(N), V(N), A(N); + U.setZero(); + V.setZero(); + A.setZero(); + solver.setU(U); + solver.setV(V); + solver.setA(A); + //! [Initialize solution] + + //! [Solve nonlinear problem] + index_t step = 50; + gsParaviewCollection collection("Deformation"); + gsMultiPatch<> displ, def; + + real_t time = 0; + real_t dt = 1e-2; + solver.setTimeStep(dt); + for (index_t k=0; k solField(def, displ); + std::string outputName = "Deformation" + std::to_string(k) + "_"; + gsWriteParaview<>( solField, outputName, 1000, true); + collection.addPart(outputName + "0.vts",time); + } + time = solver.time(); + } + //! [Solve nonlinear problem] + + // ![Save the paraview collection] + if (plot) + collection.save(); + // ![Save the paraview collection] + + return EXIT_SUCCESS; +#else + GISMO_ERROR("The tutorial needs to be compiled with gsKLShell enabled"); + return EXIT_FAILED; +#endif +}// end main \ No newline at end of file diff --git a/tutorials/nonlinear_shell_static.cpp b/tutorials/nonlinear_shell_static.cpp new file mode 100644 index 0000000..0fddd14 --- /dev/null +++ b/tutorials/nonlinear_shell_static.cpp @@ -0,0 +1,199 @@ +/** @file @file gsStructuralAnalysis/tutorials/nonlinear_shell_arclength.cpp + + @brief Tutorial for assembling the Kirchhoff-Love shell + and solving it with a static solver + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M.Verhelst +*/ + +#include + +//! [Includes] +#ifdef gsKLShell_ENABLED +#include +#endif + +#include +#include +//! [Includes] + +using namespace gismo; + +// Choose among various shell examples, default = Thin Plate +int main(int argc, char *argv[]) +{ +#ifdef gsKLShell_ENABLED + //! [Parse command line] + bool plot = false; + index_t numRefine = 1; + index_t numElevate = 1; + + gsCmdLine cmd("Nonlinear static shell tutorial."); + cmd.addInt( "e", "degreeElevation","Number of degree elevation steps to perform before solving (0: equalize degree in all directions)", numElevate ); + cmd.addInt( "r", "uniformRefine", "Number of Uniform h-refinement steps to perform before solving", numRefine ); + cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + //! [Parse command line] + + //! [Read geometry] + // Initialize [ori]ginal and [def]ormed geometry + gsMultiPatch<> ori; + gsReadFile<>("surfaces/paraboloid.xml", ori); + ori.computeTopology(); + //! [Read geometry] + + //! [Initialize geometry] + // p-refine + if (numElevate!=0) + ori.degreeElevate(numElevate); + // h-refine + for (int r =0; r < numRefine; ++r) + ori.uniformRefine(); + //! [Initialize geometry] + + //! [Construct basis] + gsMultiBasis<> bases(ori); + //! [Construct basis] + + //! [Set boundary conditions] + // Define the boundary conditions object + gsBoundaryConditions<> bc; + // Set the geometry map for computation of the Dirichlet BCs + bc.setGeoMap(ori); + + // Set the boundary conditions + bc.addCornerValue(boundary::southwest, 0.0, 0, -1); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::southeast, 0.0, 0, -1); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northwest, 0.0, 0, -1); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northeast, 0.0, 0, -1); // (corner,value, patch, unknown) + + gsPointLoads pLoads = gsPointLoads(); + gsVector<> point(2); + gsVector<> load (3); + point<< 0.5, 0.5 ; load << 0, 0.0, -1e4 ; + pLoads.addLoad(point, load, 0 ); + //! [Set boundary conditions] + + //! [Set surface force] + // The surface force is defined in the physical space, i.e. 3D + gsFunctionExpr<> force("0","0","0",3); + //! [Set surface force] + + + //! [Define the material matrix class] + // Define material parameters + // The material parameters are defined in the physical domain as well (!) + gsConstantFunction<> E (1.E9,3); + gsConstantFunction<> nu(0.45,3); + gsConstantFunction<> t (1E-2,3); + + // Define a linear material, see \ref gsMaterialMatrixLinear.h + // The first parameter is the physical domain dimension + + // gsMaterialMatrixLinear<3,real_t> materialMatrix(ori,t,E,nu); + + std::vector*> parameters{&E,&nu}; + gsOptionList options; + options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",1); + options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",1); + gsMaterialMatrixBase * materialMatrix = getMaterialMatrix<3,real_t>(ori,t,parameters,options); + //! [Define the material matrix class] + + //! [Define the assembler] + // Define the assembler. + // The first parameter is the physical domain dimension + // The second parameter is the real type + // The third parameter controls the bending stiffness + gsThinShellAssembler<3, real_t, true > assembler(ori,bases,bc,force,materialMatrix); + assembler.setPointLoads(pLoads); + //! [Define the assembler] + + + //! [Define nonlinear residual functions] + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler](gsVector const &x, gsSparseMatrix & m) + { + gsMultiPatch<> def; + assembler.constructSolution(x,def); + ThinShellAssemblerStatus status = assembler.assembleMatrix(def); + m = assembler.matrix(); + return status==ThinShellAssemblerStatus::Success; + }; + // Function for the Residual + gsStructuralAnalysisOps::Residual_t Residual = [&assembler](gsVector const &x, gsVector & v) + { + gsMultiPatch<> def; + assembler.constructSolution(x,def); + ThinShellAssemblerStatus status = assembler.assembleVector(def); + v = assembler.rhs(); + return status==ThinShellAssemblerStatus::Success; + }; + //! [Define nonlinear residual functions] + + + //! [Assemble linear part] + assembler.assemble(); + gsSparseMatrix<> K = assembler.matrix(); + gsVector<> F = assembler.rhs(); + //! [Assemble linear part] + + //! [Set static solver] + gsStaticNewton solver(K,F,Jacobian,Residual); + solver.options().setInt("verbose",2); + solver.initialize(); + //! [Set static solver] + + //! [Solve nonlinear problem] + gsInfo<<"Solving system with "< solVector = solver.solution(); + //! [Solve nonlinear problem] + + //! [Construct solution and deformed geometry] + gsMultiPatch<> displ = assembler.constructDisplacement(solVector); + gsMultiPatch<> def = assembler.constructSolution(solVector); + //! [Construct solution and deformed geometry] + + //! [Evaluate solution] + // Evaluate the solution on a reference point (parametric coordinates) + // The reference points are stored column-wise + gsMatrix<> refPoint(2,1); + refPoint<<0.5,0.5; + // Compute the values + gsVector<> physpoint = def.patch(0).eval(refPoint); + gsVector<> refVals = displ.patch(0).eval(refPoint); + // gsInfo << "Displacement at reference point: "< solField(def, displ); + gsInfo<<"Plotting in Paraview...\n"; + gsWriteParaview<>( solField, "Deformation", 1000, true); + + // Plot the membrane Von-Mises stresses on the geometry + gsPiecewiseFunction<> VMm; + assembler.constructStress(def,VMm,stress_type::von_mises_membrane); + gsField<> membraneStress(def,VMm, true); + gsWriteParaview(membraneStress,"MembraneStress",1000); + } + // ! [Export visualization in ParaView] + + return EXIT_SUCCESS; +#else + GISMO_ERROR("The tutorial needs to be compiled with gsKLShell enabled"); + return EXIT_FAILED; +#endif +}// end main