From 5a633516f6c788fd674c53372c87f5fa3428b751 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Mon, 16 Sep 2024 15:29:01 +0200 Subject: [PATCH 1/9] add static solver using optimizer --- examples/example_DisplacementControl.cpp | 378 ++++++++++++++++++ examples/example_StaticSolvers.cpp | 65 ++- src/gsStaticSolvers/gsStaticBase_.cpp | 15 + src/gsStaticSolvers/gsStaticDR.hpp | 1 - src/gsStaticSolvers/gsStaticNewton.h | 2 - src/gsStaticSolvers/gsStaticNewton.hpp | 1 - src/gsStaticSolvers/gsStaticOpt.h | 243 +++++------ src/gsStaticSolvers/gsStaticOpt.hpp | 327 ++++++++------- .../gsStructuralAnalysisTypes.h | 3 + 9 files changed, 748 insertions(+), 287 deletions(-) create mode 100644 examples/example_DisplacementControl.cpp diff --git a/examples/example_DisplacementControl.cpp b/examples/example_DisplacementControl.cpp new file mode 100644 index 0000000..b88351e --- /dev/null +++ b/examples/example_DisplacementControl.cpp @@ -0,0 +1,378 @@ +/** @file example_StaticSolvers.cpp + + @brief Demonstrations of static solvers on a 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 (2019-..., TU Delft) +*/ + +#include + +#ifdef gsKLShell_ENABLED +#include +#include +#endif + +#include +#include +#include +#include +#include + +using namespace gismo; + +void writeToCSVfile(std::string name, gsMatrix<> matrix) +{ + std::ofstream file(name.c_str()); + for(int i = 0; i < matrix.rows(); i++){ + for(int j = 0; j < matrix.cols(); j++){ + std::string str = std::to_string(matrix(i,j)); + if(j+1 == matrix.cols()){ + file< fd(assemberOptionsFile); + gsOptionList opts; + fd.getFirst(opts); + + gsMultiPatch<> mp; + gsMultiPatch<> mp_def; + + E_modulus = 1; + // thickness = 0.15; + thickness = 1; + if (!Compressibility) + PoissonRatio = 0.499; + else + PoissonRatio = 0.45; + + E_modulus = 1; + real_t bDim = 1; + real_t aDim = 1; + + mp.addPatch( gsNurbsCreator<>::BSplineSquare(1) ); // degree + mp.patch(0).coefs().col(0) *= aDim; + mp.patch(0).coefs().col(1) *= bDim; + mp.addAutoBoundaries(); + mp.embed(3); + + // p-refine + if (numElevate!=0) + mp.degreeElevate(numElevate); + + // h-refine + for (int r =0; r < numRefine; ++r) + mp.uniformRefine(); + + mp_def = mp; + gsWriteParaview<>( mp_def , "mp", 1000, true); + + //! [Refinement] + gsMultiBasis<> dbasis(mp); + + gsInfo << "Patches: "<< mp.nPatches() <<", degree: "<< dbasis.minCwiseDegree() <<"\n"; + gsInfo< bc; + bc.setGeoMap(mp); + + gsConstantFunction<> displ(0.0,3); + + gsPointLoads pLoads = gsPointLoads(); + bc.addCondition(boundary::west, condition_type::dirichlet, 0, 0, false, 0 ); // unknown 2 - z + bc.addCondition(boundary::west, condition_type::dirichlet, 0, 0, false, 1 ); // unknown 2 - z + bc.addCondition(boundary::west, condition_type::dirichlet, 0, 0, false, 2 ); // unknown 2 - z + + // bc.addCondition(boundary::north, condition_type::dirichlet, 0, 0, false, 2 ); // unknown 2 - z + // bc.addCondition(boundary::south, condition_type::dirichlet, 0, 0, false, 2 ); // unknown 2 - z + bc.addCondition(boundary::east, condition_type::dirichlet, 0, 0 ,false,1); + + bc.addCondition(boundary::east, condition_type::dirichlet, 0, 0 ,false,2); + bc.addCondition(boundary::east, condition_type::dirichlet, &displ, 0 ,false,0); + + gsVector<> tmp(3); + tmp<<0,0,0; + + //! [Refinement] + 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); + gsConstantFunction<> ratio(Ratio,3); + + real_t mu = E_modulus / (2 * (1 + PoissonRatio)); + gsConstantFunction<> alpha1(1.3,3); + gsConstantFunction<> mu1(6.3e5/4.225e5*mu,3); + gsConstantFunction<> alpha2(5.0,3); + gsConstantFunction<> mu2(0.012e5/4.225e5*mu,3); + gsConstantFunction<> alpha3(-2.0,3); + gsConstantFunction<> mu3(-0.1e5/4.225e5*mu,3); + + std::vector*> parameters; + if (material==0) // SvK & Composites + { + parameters.resize(2); + parameters[0] = &E; + parameters[1] = ν + } + else if (material==1 || material==2) // NH & NH_ext + { + parameters.resize(2); + parameters[0] = &E; + parameters[1] = ν + } + else if (material==3) // MR + { + parameters.resize(3); + parameters[0] = &E; + parameters[1] = ν + parameters[2] = ∶ + } + else if (material==4) // OG + { + parameters.resize(8); + parameters[0] = &E; + parameters[1] = ν + parameters[2] = &mu1; + parameters[3] = &alpha1; + parameters[4] = &mu2; + parameters[5] = &alpha2; + parameters[6] = &mu3; + parameters[7] = &alpha3; + } + + gsMaterialMatrixBase* materialMatrix; + + gsOptionList options; + if (material==0 && impl==1) + { + parameters.resize(2); + 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); + } + else + { + options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",material); + options.addSwitch("Compressibility","Compressibility: (false): Imcompressible | (true): Compressible",Compressibility); + options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",impl); + } + materialMatrix = getMaterialMatrix<3,real_t>(mp,t,parameters,rho,options); + + gsThinShellAssemblerBase* assembler; + assembler = new gsThinShellAssembler<3, real_t, true >(mp,dbasis,bc,force,materialMatrix); + + opts.addReal("WeakDirichlet","Penalty parameter weak dirichlet conditions",1e5); + opts.addReal("WeakClamped","Penalty parameter weak clamped conditions",1e5); + // Construct assembler object + assembler->setOptions(opts); + assembler->setPointLoads(pLoads); + + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&mp_def](gsVector const &x, gsSparseMatrix & m) + { + ThinShellAssemblerStatus status; + assembler->constructSolution(x,mp_def); + status = assembler->assembleMatrix(mp_def); + m = assembler->matrix(); + return status == ThinShellAssemblerStatus::Success; + }; + + gsStructuralAnalysisOps::Residual_t Residual = [&assembler,&mp_def](gsVector const &x, gsVector & result) + { + ThinShellAssemblerStatus status; + assembler->constructSolution(x,mp_def); + status = assembler->assembleVector(mp_def); + result = assembler->rhs(); + return status == ThinShellAssemblerStatus::Success; + }; + + gsStructuralAnalysisOps::ALResidual_t ALResidual = [&displ,&bc,&assembler,&mp_def](gsVector const &x, real_t lam, gsVector & result) + { + ThinShellAssemblerStatus status; + displ.setValue(lam,3); + assembler->updateBCs(bc); + assembler->constructSolution(x,mp_def); + status = assembler->assembleVector(mp_def); + result = assembler->rhs(); + return status == ThinShellAssemblerStatus::Success; + }; + + displ.setValue(1.0,3); + assembler->updateBCs(bc); + + assembler->assemble(); + gsSparseMatrix<> K = assembler->matrix(); + gsVector<> F = assembler->rhs(); + assembler->assembleMass(true); + gsVector<> M = assembler->rhs(); + + + gsStaticDR DRM(M,F,Residual); + gsOptionList DROptions = DRM.options(); + DROptions.setReal("damping",damping); + DROptions.setReal("alpha",alpha); + DROptions.setInt("maxIt",maxIt); + DROptions.setReal("tol",1e-2); + DROptions.setReal("tolE",1e-4); + DROptions.setInt("verbose",verbose); + DRM.setOptions(DROptions); + DRM.initialize(); + + + gsStaticNewton NWT(K,F,Jacobian,Residual); + gsOptionList NWTOptions = NWT.options(); + NWTOptions.setInt("maxIt",maxIt); + NWTOptions.setReal("tol",1e-6); + NWTOptions.setInt("verbose",verbose); + NWT.setOptions(NWTOptions); + NWT.initialize(); + + + gsStaticOpt OPT(Residual,assembler->numDofs()); + gsOptionList OPTOptions = OPT.options(); + OPTOptions.setInt("maxIt",maxIt); + OPTOptions.setReal("tol",1e-6); + OPTOptions.setInt("verbose",verbose); + OPT.setOptions(OPTOptions); + OPT.initialize(); + + gsControlDisplacement controlDR(&DRM); + gsControlDisplacement controlDC(&NWT); + gsControlDisplacement controlOPT(&OPT); + + gsStaticComposite DRNWT({&DRM,&NWT}); + DRNWT.initialize(); + DRNWT.solve(); + gsDebugVar(DRNWT.solution().norm()); + + gsStaticComposite DROPT({&DRM,&OPT}); + DROPT.initialize(); + DROPT.solve(); + gsDebugVar(DROPT.solution().norm()); + + return 0; + + // control.step(0.5); + // control.step(0.5); + controlDR.step(1.0); + // controlOPT.step(1.0); + + gsDebugVar(controlDR.solutionU().norm()); + + controlDC.step(1.0); + gsDebugVar(controlDC.solutionU().norm()); + + OPT.setUpdate(DRM.update()); + controlOPT.step(1.0); + gsDebugVar(controlOPT.solutionU().norm()); + + NWT.reset(); + NWT.setUpdate(DRM.update()); + controlDC.step(1.0); + gsDebugVar(controlDC.solutionU().norm()); + + gsVector<> displacements = controlDR.solutionU(); + + + + mp_def = assembler->constructSolution(displacements); + gsMultiPatch<> deformation = mp_def; + for (size_t k = 0; k != mp_def.nPatches(); ++k) + deformation.patch(k).coefs() -= mp.patch(k).coefs(); + + // ! [Export visualization in ParaView] + if (plot) + { + gsField<> solField(mp_def, deformation); + gsInfo<<"Plotting in Paraview...\n"; + gsWriteParaview<>( solField, "solution", 1000, true); + // ev.options().setSwitch("plot.elements", true); + // ev.writeParaview( u_sol , G, "solution"); + + // gsFileManager::open("solution.pvd"); + + gsInfo <<"Maximum deformation coef: " + << deformation.patch(0).coefs().colwise().maxCoeff() <<".\n"; + gsInfo <<"Minimum deformation coef: " + << deformation.patch(0).coefs().colwise().minCoeff() <<".\n"; + } + + delete assembler; + delete materialMatrix; + + return EXIT_SUCCESS; + +}// end main +#else//gsKLShell_ENABLED +int main(int argc, char *argv[]) +{ + gsWarn<<"G+Smo is not compiled with the gsKLShell module."; + return EXIT_FAILURE; +} +#endif diff --git a/examples/example_StaticSolvers.cpp b/examples/example_StaticSolvers.cpp index 55c4c12..089bcea 100644 --- a/examples/example_StaticSolvers.cpp +++ b/examples/example_StaticSolvers.cpp @@ -20,7 +20,10 @@ #include #include +#include +#include #include + using namespace gismo; void writeToCSVfile(std::string name, gsMatrix<> matrix) @@ -240,15 +243,13 @@ int main(int argc, char *argv[]) return status == ThinShellAssemblerStatus::Success; }; - gsStructuralAnalysisOps::ALResidual_t ALResidual = [&displ,&bc,&assembler,&mp_def](gsVector const &x, real_t lam, gsVector & result) + gsStructuralAnalysisOps::Residual_t Residual = [&assembler,&mp_def](gsVector const &x, gsVector & result) { - ThinShellAssemblerStatus status; - displ.setValue(lam,3); - assembler->updateBCs(bc); - assembler->constructSolution(x,mp_def); - status = assembler->assembleVector(mp_def); - result = assembler->rhs(); - return status == ThinShellAssemblerStatus::Success; + ThinShellAssemblerStatus status; + assembler->constructSolution(x,mp_def); + status = assembler->assembleVector(mp_def); + result = assembler->rhs(); + return status == ThinShellAssemblerStatus::Success; }; displ.setValue(1.0,3); @@ -261,7 +262,7 @@ int main(int argc, char *argv[]) gsVector<> M = assembler->rhs(); - gsStaticDR DRM(M,F,ALResidual); + gsStaticDR DRM(M,F,Residual); gsOptionList DROptions = DRM.options(); DROptions.setReal("damping",damping); DROptions.setReal("alpha",alpha); @@ -270,28 +271,52 @@ int main(int argc, char *argv[]) DROptions.setReal("tolE",1e-4); DROptions.setInt("verbose",verbose); DRM.setOptions(DROptions); + DRM.initialize(); - gsStaticNewton NWT(K,F,Jacobian,ALResidual); + gsStaticNewton NWT(K,F,Jacobian,Residual); gsOptionList NWTOptions = NWT.options(); NWTOptions.setInt("maxIt",maxIt); NWTOptions.setReal("tol",1e-6); NWTOptions.setInt("verbose",verbose); NWT.setOptions(NWTOptions); + NWT.initialize(); + + + gsStaticOpt OPT(Residual,assembler->numDofs()); + gsOptionList OPTOptions = OPT.options(); + OPTOptions.setInt("maxIt",maxIt); + OPTOptions.setReal("tol",1e-6); + OPTOptions.setInt("verbose",verbose); + OPT.setOptions(OPTOptions); + OPT.initialize(); + + // gsInfo<<"Dynamic relaxation\n"; + // DRM.solve(); + // gsInfo<<"Solution: "< controlDR(&DRM); - gsControlDisplacement controlDC(&NWT); - // control.step(0.5); - // control.step(0.5); - controlDR.step(1.0); + gsStaticComposite DRNWT({&DRM,&NWT}); + DRNWT.initialize(); + gsInfo<<"Dynamic Relaxation -> Newton\n"; + DRNWT.solve(); + gsInfo<<"Solution: "< DROPT({&DRM,&OPT}); + DROPT.initialize(); + gsInfo<<"Dynamic Relaxation -> Optimizer\n"; + DROPT.solve(); + gsInfo<<"Solution: "< displacements = controlDR.solutionU(); + gsVector<> displacements = NWT.solution(); mp_def = assembler->constructSolution(displacements); gsMultiPatch<> deformation = mp_def; diff --git a/src/gsStaticSolvers/gsStaticBase_.cpp b/src/gsStaticSolvers/gsStaticBase_.cpp index c00e201..d153fe3 100644 --- a/src/gsStaticSolvers/gsStaticBase_.cpp +++ b/src/gsStaticSolvers/gsStaticBase_.cpp @@ -9,6 +9,14 @@ #include #include +#include +#include + +#include +#ifdef gsHLBFGS_ENABLED +#include +#endif + #include #include @@ -18,4 +26,11 @@ namespace gismo CLASS_TEMPLATE_INST gsStaticDR; CLASS_TEMPLATE_INST gsStaticNewton; CLASS_TEMPLATE_INST gsStaticComposite; + + CLASS_TEMPLATE_INST gsOptProblemStatic; + CLASS_TEMPLATE_INST gsStaticOpt>; +#ifdef gsHLBFGS_ENABLED + CLASS_TEMPLATE_INST gsStaticOpt>; +#endif + } diff --git a/src/gsStaticSolvers/gsStaticDR.hpp b/src/gsStaticSolvers/gsStaticDR.hpp index 6a4dc6b..10240a9 100644 --- a/src/gsStaticSolvers/gsStaticDR.hpp +++ b/src/gsStaticSolvers/gsStaticDR.hpp @@ -11,7 +11,6 @@ Author(s): H.M. Verhelst (2019-..., TU Delft) */ -#include #pragma once namespace gismo diff --git a/src/gsStaticSolvers/gsStaticNewton.h b/src/gsStaticSolvers/gsStaticNewton.h index 80ad71f..6b2b8ff 100644 --- a/src/gsStaticSolvers/gsStaticNewton.h +++ b/src/gsStaticSolvers/gsStaticNewton.h @@ -11,8 +11,6 @@ Author(s): H.M. Verhelst (2019-..., TU Delft) */ -#include - #include #pragma once diff --git a/src/gsStaticSolvers/gsStaticNewton.hpp b/src/gsStaticSolvers/gsStaticNewton.hpp index 91921ef..0a47c6a 100644 --- a/src/gsStaticSolvers/gsStaticNewton.hpp +++ b/src/gsStaticSolvers/gsStaticNewton.hpp @@ -11,7 +11,6 @@ Author(s): H.M. Verhelst (2019-..., TU Delft) */ -#include #pragma once namespace gismo diff --git a/src/gsStaticSolvers/gsStaticOpt.h b/src/gsStaticSolvers/gsStaticOpt.h index faf0123..7800905 100644 --- a/src/gsStaticSolvers/gsStaticOpt.h +++ b/src/gsStaticSolvers/gsStaticOpt.h @@ -12,10 +12,91 @@ */ #include +#include + +#pragma once namespace gismo { +template +class gsOptProblemStatic : public gsOptProblem +{ +protected: + + typedef gsOptProblem Base; + typedef typename gsStructuralAnalysisOps::Residual_t Residual_t; + typedef typename gsStructuralAnalysisOps::ALResidual_t ALResidual_t; + +public: + + gsOptProblemStatic(const typename gsStructuralAnalysisOps::Residual_t & residualFun, const T & L, index_t numDesignVars) + : + m_residualFun(residualFun), + m_L(L) // unused, in fact, but needs to be assigned + { + m_numDesignVars = numDesignVars; + m_curDesign.resize(numDesignVars,1); + m_curDesign.setZero(); + } + + gsOptProblemStatic(const typename gsStructuralAnalysisOps::ALResidual_t & ALresidualFun, const T & L, index_t numDesignVars) + : + m_ALresidualFun(ALresidualFun), + m_L(L) + { + m_numDesignVars = numDesignVars; + m_residualFun = [this](gsVector const & x, gsVector & result) -> bool + { + return m_ALresidualFun(x,m_L,result); + }; + m_curDesign.resize(numDesignVars,1); + m_curDesign.setZero(); + } + +public: + + T evalObj( const gsAsConstVector & u ) const + { + gsVector result; + m_residualFun(u,result); + return result.norm(); + } + + void gradObj_into( const gsAsConstVector & u, gsAsVector & result) const + { + gsVector tmp; + result.resize(u.rows()); + m_residualFun(u,tmp); + result = -tmp; + } + + // void evalCon_into( const gsAsConstVector & u, gsAsVector & result) const; + + // void jacobCon_into( const gsAsConstVector & u, gsAsVector & result) const; + +private: + + Residual_t m_residualFun; + ALResidual_t m_ALresidualFun; + const T & m_L; + // Lastly, we forward the memebers of the base clase gsOptProblem + using Base::m_numDesignVars; + using Base::m_numConstraints; + using Base::m_numConJacNonZero; + + using Base::m_desLowerBounds; + using Base::m_desUpperBounds; + + using Base::m_conLowerBounds; + using Base::m_conUpperBounds; + + using Base::m_conJacRows; + using Base::m_conJacCols; + + using Base::m_curDesign; +}; + /** * @brief Static solver using the Dynamic Relaxation method * @@ -23,7 +104,7 @@ namespace gismo * * \ingroup gsStaticBase */ -template +template > class gsStaticOpt : public gsStaticBase { protected: @@ -40,137 +121,109 @@ class gsStaticOpt : public gsStaticBase * * @param[in] Residual The residual */ - gsStaticOpt( const Residual_t &Residual ) + gsStaticOpt( const Residual_t &Residual, const index_t numDofs ) : - m_residualFun(Residual), - m_ALresidualFun(nullptr), - m_jacobian(nullptr) + m_optProblem(Residual,m_L,numDofs), + m_optimizer(&m_optProblem) { + m_dofs = numDofs; this->_init(); } - /** - * @brief Constructor - * - * @param[in] Residual The residual - * @param[in] Jacobian The jacobian - */ - gsStaticOpt( const Residual_t &Residual, - const Jacobian_t &Jacobian ) - : - m_residualFun(Residual), - m_ALresidualFun(nullptr), - m_jacobian(Jacobian) - { - this->_init(); - } - - /** * @brief Constructs a new instance. * * @param[in] ALResidual The residual as arc-length object - * @param[in] Jacobian The jacobian */ - gsStaticOpt( const ALResidual_t & ALResidual ) + gsStaticOpt( const ALResidual_t & ALResidual, const index_t numDofs ) : - m_ALresidualFun(ALResidual), - m_jacobian(nullptr) + m_optProblem(ALResidual,m_L,numDofs), + m_optimizer(&m_optProblem) { m_L = 1.0; - m_residualFun = [this](gsVector const & x, gsVector & result) -> bool - { - return m_ALresidualFun(x,m_L,result); - }; + m_dofs = numDofs; this->_init(); } - /** - * @brief Constructs a new instance. - * - * @param[in] ALResidual The residual as arc-length object - * @param[in] Jacobian The jacobian - */ - gsStaticOpt( const ALResidual_t & ALResidual, - const Jacobian_t &Jacobian ) - : - m_ALresidualFun(ALResidual), - m_jacobian(Jacobian) - { - m_L = 1.0; - m_residualFun = [this](gsVector const & x, gsVector & result) -> bool - { - return m_ALresidualFun(x,m_L,result); - }; - this->_init(); - } +public: + gsOptionList & optimizerOptions() { return m_optimizer.options(); } /// gsStaticBase base functions public: /// See \ref gsStaticBase gsStatus solve() override; - /// See \ref gsStaticBase - void initialize() override; + // /// See \ref gsStaticBase + // void initialize() override; - /// See \ref gsStaticBase - void initOutput() override; + // /// See \ref gsStaticBase + // void initOutput() override; - /// See \ref gsStaticBase - void stepOutput(index_t k) override; + // /// See \ref gsStaticBase + // void stepOutput(index_t k) override; /// See \ref gsStaticBase void defaultOptions() override; - /// See \ref gsStaticBase - void reset() override; + // /// See \ref gsStaticBase + // void reset() override; /// See \ref gsStaticBase void getOptions() override; + /// Returns the stability indicator + T indicator(const gsSparseMatrix &, T) + { + GISMO_NO_IMPLEMENTATION; + } + + /// Returns the stability vector + gsVector stabilityVec(const gsSparseMatrix &, T) + { + GISMO_NO_IMPLEMENTATION; + } + +// Class-specific functions +// public: +// void initialize(); + protected: /// See \ref solve() void _solve(); /// Initializes the method void _init(); - gsVector _computeResidual(const gsVector & U); +// gsVector _computeResidual(const gsVector & U); -public: +// public: - /// Return the residual norm - T residualNorm() const { return m_R.norm(); } +// /// Return the residual norm +// T residualNorm() const { return m_R.norm(); } protected: - Residual_t m_residualFun; - const ALResidual_t m_ALresidualFun; - const Jacobian_t m_jacobian; + Optimizer m_optimizer; + gsOptProblemStatic m_optProblem; - // Solution + // // Solution using Base::m_U; using Base::m_DeltaU; - using Base::m_deltaU; + // using Base::m_deltaU; using Base::m_L; - using Base::m_DeltaL; - using Base::m_deltaL; + // using Base::m_DeltaL; + // using Base::m_deltaL; // Iterations using Base::m_numIterations; using Base::m_maxIterations; - // Residuals - using Base::m_R; + // // Residuals + // using Base::m_R; - // Tolerances + // // Tolerances using Base::m_tolF; using Base::m_tolU; - // Residual norms - using Base::m_residual; - using Base::m_residualIni; - using Base::m_residualOld; - // Options using Base::m_options; using Base::m_verbose; @@ -185,48 +238,10 @@ class gsStaticOpt : public gsStaticBase using Base::m_status; }; -//! [OptProblemExample Class] -template -class gsOptProblemStatic : public gsOptProblem -//! [OptProblemExample Class] -{ -public: - - gsOptProblemStatic(); - -public: - - T evalObj( const gsAsConstVector & u ) const; - void gradObj_into( const gsAsConstVector & u, gsAsVector & result) const; - - void evalCon_into( const gsAsConstVector & u, gsAsVector & result) const; - - void jacobCon_into( const gsAsConstVector & u, gsAsVector & result) const; - -private: - - // Lastly, we forward the memebers of the base clase gsOptProblem - using gsOptProblem::m_numDesignVars; - using gsOptProblem::m_numConstraints; - using gsOptProblem::m_numConJacNonZero; - - using gsOptProblem::m_desLowerBounds; - using gsOptProblem::m_desUpperBounds; - - using gsOptProblem::m_conLowerBounds; - using gsOptProblem::m_conUpperBounds; - - using gsOptProblem::m_conJacRows; - using gsOptProblem::m_conJacCols; - - using gsOptProblem::m_curDesign; -}; -//! [OptProblem] } //namespace - #ifndef GISMO_BUILD_LIB #include GISMO_HPP_HEADER(gsStaticOpt.hpp) #endif diff --git a/src/gsStaticSolvers/gsStaticOpt.hpp b/src/gsStaticSolvers/gsStaticOpt.hpp index 280bb43..943a3c2 100644 --- a/src/gsStaticSolvers/gsStaticOpt.hpp +++ b/src/gsStaticSolvers/gsStaticOpt.hpp @@ -11,7 +11,6 @@ Author(s): H.M. Verhelst (2019-..., TU Delft) */ -#include #include #ifdef gsHLBFGS_ENABLED @@ -27,176 +26,206 @@ namespace gismo { -template -void gsStaticOpt::defaultOptions() +template +void gsStaticOpt::defaultOptions() { Base::defaultOptions(); } -template -void gsStaticOpt::getOptions() +template +void gsStaticOpt::_init() { - Base::getOptions(); -} + m_U.setZero(m_dofs); + m_DeltaU.setZero(m_dofs); -template -void gsStaticOpt::initOutput() -{ - gsInfo<<"\t"; - gsInfo< -void gsStaticOpt::stepOutput(index_t k) -{ - gsInfo<<"\t"; - gsInfo<(m_residualFun, m_dofs); + // m_optimizer = Optimizer(&m_optProblem); } -template -gsStatus gsStaticOpt::solve() +template +void gsStaticOpt::getOptions() { - try - { - _solve(); - m_status = gsStatus::Success; - } - catch (int errorCode) + Base::getOptions(); + + m_optimizer.options().setInt("MaxIterations",m_maxIterations); + m_optimizer.options().setInt("Verbose",m_verbose); + if (dynamic_cast*>(&m_optimizer)) { - 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; + m_optimizer.options().setReal("MinGradientLength",m_tolF); + m_optimizer.options().setReal("MinStepLength",m_tolU); } - catch (...) +#ifdef gsHLBFGS_ENABLED + else if (dynamic_cast*>(&m_optimizer)) { - m_status = gsStatus::OtherError; + m_optimizer.options().setReal("MinGradientLength",m_tolF); + m_optimizer.options().setReal("MinStepLength",m_tolU); } - return m_status; -} - -template -void gsStaticOpt::_solve() -{ - m_Eks.clear(); - m_Eks.reserve(m_maxIterations); - - if (m_verbose) initOutput(); - - _start(); - - m_Ek0 = m_Ek; - m_Eks.push_back(m_Ek); - - if (m_verbose != 0) stepOutput(0); - index_t resetIt = 0; - for (m_numIterations=1; m_numIterations!=m_maxIterations; m_numIterations++, resetIt++) +#endif +#ifdef gsIpOpt_ENABLED + else if (dynamic_cast*>(&m_optimizer)) { - _iteration(); - if ((m_c==0 && m_Ek_prev > m_Ek) || resetIt==m_resetIterations)// || (m_Ek/m_Ek_prev > 1/m_tolE && m_Ek_prev!=0)) - { - resetIt = 0; - _peak(); - } - - if (m_verbose!=0) - if (m_numIterations % m_verbose == 0 || m_verbose==-1 ) stepOutput(m_numIterations); - - m_Eks.push_back(m_Ek); - - m_residualOld = m_residual; - - if (m_residual/m_residualIni < m_tolF && m_Ek/m_Ek0 < m_tolE && m_deltaU.norm()/m_DeltaU.norm() < m_tolU) - { - m_U += m_DeltaU; - gsDebug <<"Converged: \n"; - gsDebug <<"\t |R|/|R0| = "< -void gsStaticOpt::initialize() -{ - this->reset(); - getOptions(); -} - -template -void gsStaticOpt::reset() -{ - Base::reset(); - m_status = gsStatus::NotStarted; +#endif + else + GISMO_ERROR("Optimizer not recognized\n"); } -template -void gsStaticOpt::_init() +template +gsStatus gsStaticOpt::solve() { - this->reset(); - defaultOptions(); + this->getOptions(); + m_optimizer.solve(m_U+m_DeltaU); + + m_DeltaU = m_optimizer.currentDesign() - m_U; + m_U += m_DeltaU; + + m_numIterations = m_optimizer.iterations(); + return gsStatus::Success; } -template -gsVector gsStaticOpt::_computeResidual(const gsVector & U) -{ - gsVector resVec; - if (!m_residualFun(U, resVec)) - throw 2; - return resVec; -} -template -gsOptProblemStatic::gsOptProblemStatic() -{ - // Number of design variables - m_numDesignVars = 2; - // Number of constraints - m_numConstraints = 0; - // Number of non-zeros in the Jacobian of the constraints - m_numConJacNonZero = 0; - - m_desLowerBounds.resize(m_numDesignVars); - m_desUpperBounds.resize(m_numDesignVars); - - m_desLowerBounds.setConstant(-1e19); - m_desUpperBounds.setConstant( 1e19); - - // we initialize x in bounds, in the upper right quadrant - m_curDesign.resize(m_numDesignVars,1); - m_curDesign.setZero(); -} +//! [OptProblem] + +// template +// void gsStaticOpt::getOptions() +// { +// Base::getOptions(); +// } + +// template +// gsStatus gsStaticOpt::solve() +// { +// try +// { +// _solve(); +// 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 gsStaticOpt::_solve() +// { +// m_Eks.clear(); +// m_Eks.reserve(m_maxIterations); + +// if (m_verbose) initOutput(); + +// _start(); + +// m_Ek0 = m_Ek; +// m_Eks.push_back(m_Ek); + +// if (m_verbose != 0) stepOutput(0); +// index_t resetIt = 0; +// for (m_numIterations=1; m_numIterations!=m_maxIterations; m_numIterations++, resetIt++) +// { +// _iteration(); +// if ((m_c==0 && m_Ek_prev > m_Ek) || resetIt==m_resetIterations)// || (m_Ek/m_Ek_prev > 1/m_tolE && m_Ek_prev!=0)) +// { +// resetIt = 0; +// _peak(); +// } + +// if (m_verbose!=0) +// if (m_numIterations % m_verbose == 0 || m_verbose==-1 ) stepOutput(m_numIterations); + +// m_Eks.push_back(m_Ek); + +// m_residualOld = m_residual; + +// if (m_residual/m_residualIni < m_tolF && m_Ek/m_Ek0 < m_tolE && m_deltaU.norm()/m_DeltaU.norm() < m_tolU) +// { +// m_U += m_DeltaU; +// gsDebug <<"Converged: \n"; +// gsDebug <<"\t |R|/|R0| = "< +// void gsStaticOpt::initialize() +// { +// this->reset(); +// getOptions(); +// } + +// template +// void gsStaticOpt::reset() +// { +// Base::reset(); +// m_status = gsStatus::NotStarted; +// } + +// template +// void gsStaticOpt::_init() +// { +// this->reset(); +// defaultOptions(); +// } + +// template +// gsVector gsStaticOpt::_computeResidual(const gsVector & U) +// { +// gsVector resVec; +// if (!m_residualFun(U, resVec)) +// throw 2; +// return resVec; +// } + +// template +// gsOptProblemStatic::gsOptProblemStatic() +// { + +// // Number of design variables +// m_numDesignVars = 2; +// // Number of constraints +// m_numConstraints = 0; +// // Number of non-zeros in the Jacobian of the constraints +// m_numConJacNonZero = 0; + +// m_desLowerBounds.resize(m_numDesignVars); +// m_desUpperBounds.resize(m_numDesignVars); + +// m_desLowerBounds.setConstant(-1e19); +// m_desUpperBounds.setConstant( 1e19); + +// // we initialize x in bounds, in the upper right quadrant +// m_curDesign.resize(m_numDesignVars,1); +// m_curDesign.setZero(); +// } } // namespace gismo diff --git a/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h b/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h index f571432..e58e8c1 100644 --- a/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h +++ b/src/gsStructuralAnalysisTools/gsStructuralAnalysisTypes.h @@ -56,6 +56,9 @@ enum struct gsStatus template struct gsStructuralAnalysisOps { + /// Residual energy + typedef std::function < bool ( gsVector const &, T &)> Energy_t; + /// Force typedef std::function < bool ( gsVector & )> Force_t; /// Time-dependent force From 5efff22826b212f4943a88faf721a1c82fd4cf73 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Mon, 16 Sep 2024 15:38:40 +0200 Subject: [PATCH 2/9] add displacement-control example --- examples/example_DisplacementControl.cpp | 77 ++++++--------------- src/gsStaticSolvers/gsControlDisplacement.h | 8 +++ 2 files changed, 29 insertions(+), 56 deletions(-) diff --git a/examples/example_DisplacementControl.cpp b/examples/example_DisplacementControl.cpp index b88351e..5643838 100644 --- a/examples/example_DisplacementControl.cpp +++ b/examples/example_DisplacementControl.cpp @@ -20,7 +20,6 @@ #include #include -#include #include #include @@ -243,15 +242,6 @@ int main(int argc, char *argv[]) return status == ThinShellAssemblerStatus::Success; }; - gsStructuralAnalysisOps::Residual_t Residual = [&assembler,&mp_def](gsVector const &x, gsVector & result) - { - ThinShellAssemblerStatus status; - assembler->constructSolution(x,mp_def); - status = assembler->assembleVector(mp_def); - result = assembler->rhs(); - return status == ThinShellAssemblerStatus::Success; - }; - gsStructuralAnalysisOps::ALResidual_t ALResidual = [&displ,&bc,&assembler,&mp_def](gsVector const &x, real_t lam, gsVector & result) { ThinShellAssemblerStatus status; @@ -273,7 +263,11 @@ int main(int argc, char *argv[]) gsVector<> M = assembler->rhs(); - gsStaticDR DRM(M,F,Residual); + + + + + gsStaticDR DRM(M,F,ALResidual); gsOptionList DROptions = DRM.options(); DROptions.setReal("damping",damping); DROptions.setReal("alpha",alpha); @@ -284,8 +278,13 @@ int main(int argc, char *argv[]) DRM.setOptions(DROptions); DRM.initialize(); + // gsControlDisplacement controlDR(&DRM); + // controlDR.step(0.5); + // controlDR.step(0.5); + // controlDR.reset(); + // controlDR.step(1.0); - gsStaticNewton NWT(K,F,Jacobian,Residual); + gsStaticNewton NWT(K,F,Jacobian,ALResidual); gsOptionList NWTOptions = NWT.options(); NWTOptions.setInt("maxIt",maxIt); NWTOptions.setReal("tol",1e-6); @@ -293,53 +292,19 @@ int main(int argc, char *argv[]) NWT.setOptions(NWTOptions); NWT.initialize(); - - gsStaticOpt OPT(Residual,assembler->numDofs()); - gsOptionList OPTOptions = OPT.options(); - OPTOptions.setInt("maxIt",maxIt); - OPTOptions.setReal("tol",1e-6); - OPTOptions.setInt("verbose",verbose); - OPT.setOptions(OPTOptions); - OPT.initialize(); - - gsControlDisplacement controlDR(&DRM); gsControlDisplacement controlDC(&NWT); - gsControlDisplacement controlOPT(&OPT); - - gsStaticComposite DRNWT({&DRM,&NWT}); - DRNWT.initialize(); - DRNWT.solve(); - gsDebugVar(DRNWT.solution().norm()); - - gsStaticComposite DROPT({&DRM,&OPT}); - DROPT.initialize(); - DROPT.solve(); - gsDebugVar(DROPT.solution().norm()); - - return 0; - - // control.step(0.5); - // control.step(0.5); - controlDR.step(1.0); - // controlOPT.step(1.0); - - gsDebugVar(controlDR.solutionU().norm()); - + gsInfo<<"Step 1 (dL=0.5)\n"; + controlDC.step(0.5); + gsInfo<<"Step 2 (dL=0.5)\n"; + controlDC.step(0.5); + gsInfo<<"Solution norm after 2 steps: "< displacements = controlDR.solutionU(); - + gsInfo<<"Solution norm after 1 step : "< displacements = controlDC.solutionU(); mp_def = assembler->constructSolution(displacements); gsMultiPatch<> deformation = mp_def; diff --git a/src/gsStaticSolvers/gsControlDisplacement.h b/src/gsStaticSolvers/gsControlDisplacement.h index 1aadea9..bec0011 100644 --- a/src/gsStaticSolvers/gsControlDisplacement.h +++ b/src/gsStaticSolvers/gsControlDisplacement.h @@ -77,6 +77,14 @@ class gsControlDisplacement : public gsContinuationBase m_solver->reset(); } + void setZero() + { + m_solver->reset(); + m_U = gsVector::Zero(m_solver->numDofs()); + m_L = 0; + m_solver->setLoad(m_L); + } + protected: mutable gsStaticBase * m_solver; T m_L; From c813be435b7c98a2e4c97be0a5396e47e81b310d Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Mon, 16 Sep 2024 18:58:12 +0200 Subject: [PATCH 3/9] add optim --- src/gsStaticSolvers/gsStaticOpt.hpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/gsStaticSolvers/gsStaticOpt.hpp b/src/gsStaticSolvers/gsStaticOpt.hpp index 943a3c2..7ca2c6a 100644 --- a/src/gsStaticSolvers/gsStaticOpt.hpp +++ b/src/gsStaticSolvers/gsStaticOpt.hpp @@ -17,6 +17,10 @@ #include #endif +#ifdef gsOptim_ENABLED +#include +#endif + #ifdef gsIpOpt_ENABLED #include #endif @@ -66,6 +70,13 @@ void gsStaticOpt::getOptions() m_optimizer.options().setReal("MinStepLength",m_tolU); } #endif +#ifdef gsOptim_ENABLED + else if (dynamic_cast*>(&m_optimizer)) + { + m_optimizer.options().setReal("RelObjFnChangeTol",m_tolF); + m_optimizer.options().setReal("RelSolChangeTol",m_tolU); + } +#endif #ifdef gsIpOpt_ENABLED else if (dynamic_cast*>(&m_optimizer)) { From 93ed51a7a2f6a89e2bbc5f3efd71f4c5daa20fe1 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 4 Dec 2024 11:10:45 +0100 Subject: [PATCH 4/9] static solver changes --- src/gsStaticSolvers/gsStaticBase_.cpp | 16 ++++++++++++++++ src/gsStaticSolvers/gsStaticOpt.h | 8 ++++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/src/gsStaticSolvers/gsStaticBase_.cpp b/src/gsStaticSolvers/gsStaticBase_.cpp index d153fe3..67759c8 100644 --- a/src/gsStaticSolvers/gsStaticBase_.cpp +++ b/src/gsStaticSolvers/gsStaticBase_.cpp @@ -16,6 +16,9 @@ #ifdef gsHLBFGS_ENABLED #include #endif +#ifdef gsOptim_ENABLED +#include +#endif #include #include @@ -32,5 +35,18 @@ namespace gismo #ifdef gsHLBFGS_ENABLED CLASS_TEMPLATE_INST gsStaticOpt>; #endif +#ifdef gsOptim_ENABLED + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + // CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; + CLASS_TEMPLATE_INST gsStaticOpt>; +#endif } diff --git a/src/gsStaticSolvers/gsStaticOpt.h b/src/gsStaticSolvers/gsStaticOpt.h index 7800905..dd224a9 100644 --- a/src/gsStaticSolvers/gsStaticOpt.h +++ b/src/gsStaticSolvers/gsStaticOpt.h @@ -123,8 +123,8 @@ class gsStaticOpt : public gsStaticBase */ gsStaticOpt( const Residual_t &Residual, const index_t numDofs ) : - m_optProblem(Residual,m_L,numDofs), - m_optimizer(&m_optProblem) + m_optimizer(&m_optProblem), + m_optProblem(Residual,m_L,numDofs) { m_dofs = numDofs; this->_init(); @@ -137,8 +137,8 @@ class gsStaticOpt : public gsStaticBase */ gsStaticOpt( const ALResidual_t & ALResidual, const index_t numDofs ) : - m_optProblem(ALResidual,m_L,numDofs), - m_optimizer(&m_optProblem) + m_optimizer(&m_optProblem), + m_optProblem(ALResidual,m_L,numDofs) { m_L = 1.0; m_dofs = numDofs; From a0dc5c333bdcbb001de1d29b55a4f91d03f3f23c Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Thu, 5 Dec 2024 16:19:39 +0100 Subject: [PATCH 5/9] add `gsStaticOpt` in example only add solution in `gsStaticComposite` if converged --- src/gsStaticSolvers/gsStaticComposite.hpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/gsStaticSolvers/gsStaticComposite.hpp b/src/gsStaticSolvers/gsStaticComposite.hpp index b84c3f6..61b7785 100644 --- a/src/gsStaticSolvers/gsStaticComposite.hpp +++ b/src/gsStaticSolvers/gsStaticComposite.hpp @@ -27,16 +27,21 @@ gsStatus gsStaticComposite::solve() { if (m_verbose > 0) gsInfo<<"Solver "<0) { m_solvers[k]->setUpdate(m_solvers[k-1]->update()); } m_status = m_solvers[k]->solve(); - m_numIterations += m_solvers[k]->iterations(); - m_U = m_solvers[k]->solution(); - m_DeltaU = m_solvers[k]->update(); + if (m_status==gsStatus::Success) + { + m_numIterations += m_solvers[k]->iterations(); + m_U = m_solvers[k]->solution(); + m_DeltaU = m_solvers[k]->update(); + } + else if (m_verbose > 0) + gsInfo<<"Solver "< Date: Sun, 8 Dec 2024 19:36:38 +0100 Subject: [PATCH 6/9] add unittest add docs static optimizer change cmakelists targets --- .github/workflows/ci.yml | 2 +- CMakeLists.txt | 8 +- benchmarks/CMakeLists.txt | 2 +- examples/CMakeLists.txt | 2 +- solvers/CMakeLists.txt | 2 +- src/gsStaticSolvers/gsStaticOpt.h | 88 ++++-- src/gsStaticSolvers/gsStaticOpt.hpp | 141 +--------- tutorials/CMakeLists.txt | 2 +- unittests/gsStaticSolver_test.cpp | 421 ++++++++++++++++++++++++++++ 9 files changed, 494 insertions(+), 174 deletions(-) create mode 100644 unittests/gsStaticSolver_test.cpp diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 21e42de..772260b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -49,4 +49,4 @@ jobs: - name: "Run for ${{ matrix.os }}" shell: bash working-directory: ${{runner.workspace}} - run: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" -D CTEST_SITE="${{ matrix.os }}_[actions]" -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON" -D GISMO_OPTIONAL='gsSpectra\\;gsKLShell\\;gsElasticity\\;${{ github.event.repository.name }}' -Q + run: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" -D CTEST_CONFIGURATION_TYPE=RelWithDebInfo -D LABELS_FOR_SUBPROJECTS="gsStructuralAnalysis-benchmarks" -D CTEST_SITE="${{ matrix.os }}_[actions]" -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON" -D GISMO_OPTIONAL='gsOptim\\;gsSpectra\\;gsKLShell\\;gsElasticity\\;${{ github.event.repository.name }}' -Q \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index cf3d4d0..19ec531 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -46,10 +46,10 @@ else() add_subdirectory(tutorials EXCLUDE_FROM_ALL) endif(GISMO_BUILD_EXAMPLES) -# # add unittests -# aux_gs_cpp_directory(${PROJECT_SOURCE_DIR}/unittests unittests_SRCS) -# set(gismo_UNITTESTS ${gismo_UNITTESTS} ${unittests_SRCS} -# CACHE INTERNAL "gismo list of unittests") +# add unittests +aux_gs_cpp_directory(${PROJECT_SOURCE_DIR}/unittests unittests_SRCS) +set(gismo_UNITTESTS ${gismo_UNITTESTS} ${unittests_SRCS} + CACHE INTERNAL "gismo list of unittests") # needed?: set(EXECUTABLE_OUTPUT_PATH ${CMAKE_BINARY_DIR}/bin/) diff --git a/benchmarks/CMakeLists.txt b/benchmarks/CMakeLists.txt index 06dcbf9..8688172 100644 --- a/benchmarks/CMakeLists.txt +++ b/benchmarks/CMakeLists.txt @@ -11,7 +11,7 @@ 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}") + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-benchmarks") if(GISMO_BUILD_EXAMPLES) set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") else(GISMO_BUILD_EXAMPLES) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 689a545..6342ca0 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -11,7 +11,7 @@ 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}") + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-examples") if(GISMO_BUILD_EXAMPLES) set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") else(GISMO_BUILD_EXAMPLES) diff --git a/solvers/CMakeLists.txt b/solvers/CMakeLists.txt index 01fe917..13be6b8 100644 --- a/solvers/CMakeLists.txt +++ b/solvers/CMakeLists.txt @@ -11,7 +11,7 @@ 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}") + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-solvers") if(GISMO_BUILD_EXAMPLES) set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") else(GISMO_BUILD_EXAMPLES) diff --git a/src/gsStaticSolvers/gsStaticOpt.h b/src/gsStaticSolvers/gsStaticOpt.h index dd224a9..6322972 100644 --- a/src/gsStaticSolvers/gsStaticOpt.h +++ b/src/gsStaticSolvers/gsStaticOpt.h @@ -15,11 +15,27 @@ #include #pragma once - namespace gismo { -template +/** + * @file gsStaticOpt.h + * @brief Header file for the gsStaticOpt and gsOptProblemStatic classes. + * + * This file contains the definitions for the gsStaticOpt and gsOptProblemStatic classes, + * which are used for static optimization in structural analysis. + */ + +/** + * @class gsOptProblemStatic + * @brief A class representing a static optimization problem. + * + * @tparam T The coefficient type. + * + * This class inherits from gsOptProblem and provides functionality for evaluating + * the objective function and its gradient for static optimization problems. + */ +template class gsOptProblemStatic : public gsOptProblem { protected: @@ -29,7 +45,13 @@ class gsOptProblemStatic : public gsOptProblem typedef typename gsStructuralAnalysisOps::ALResidual_t ALResidual_t; public: - + /** + * @brief Constructor for gsOptProblemStatic. + * + * @param residualFun The residual function. + * @param L A coefficient (unused but needs to be assigned). + * @param numDesignVars The number of design variables. + */ gsOptProblemStatic(const typename gsStructuralAnalysisOps::Residual_t & residualFun, const T & L, index_t numDesignVars) : m_residualFun(residualFun), @@ -40,6 +62,13 @@ class gsOptProblemStatic : public gsOptProblem m_curDesign.setZero(); } + /** + * @brief Constructor for gsOptProblemStatic with arc-length residual function. + * + * @param ALresidualFun The arc-length residual function. + * @param L A coefficient. + * @param numDesignVars The number of design variables. + */ gsOptProblemStatic(const typename gsStructuralAnalysisOps::ALResidual_t & ALresidualFun, const T & L, index_t numDesignVars) : m_ALresidualFun(ALresidualFun), @@ -55,7 +84,12 @@ class gsOptProblemStatic : public gsOptProblem } public: - + /** + * @brief Evaluates the objective function. + * + * @param u The input vector. + * @return The norm of the residual. + */ T evalObj( const gsAsConstVector & u ) const { gsVector result; @@ -63,6 +97,12 @@ class gsOptProblemStatic : public gsOptProblem return result.norm(); } + /** + * @brief Computes the gradient of the objective function. + * + * @param u The input vector. + * @param result The output vector for the gradient. + */ void gradObj_into( const gsAsConstVector & u, gsAsVector & result) const { gsVector tmp; @@ -71,10 +111,6 @@ class gsOptProblemStatic : public gsOptProblem result = -tmp; } - // void evalCon_into( const gsAsConstVector & u, gsAsVector & result) const; - - // void jacobCon_into( const gsAsConstVector & u, gsAsVector & result) const; - private: Residual_t m_residualFun; @@ -98,10 +134,15 @@ class gsOptProblemStatic : public gsOptProblem }; /** - * @brief Static solver using the Dynamic Relaxation method + * @class gsStaticOpt + * @brief Static solver using the Dynamic Relaxation method. * - * @tparam T coefficient type + * @tparam T The coefficient type. + * @tparam Optimizer The optimizer type (default is gsGradientDescent). * + * This class inherits from gsStaticBase and provides functionality for solving + * static optimization problems using the Dynamic Relaxation method. + * \ingroup gsStaticBase */ template > @@ -117,9 +158,10 @@ class gsStaticOpt : public gsStaticBase public: /** - * @brief Constructor + * @brief Constructor for gsStaticOpt. * - * @param[in] Residual The residual + * @param Residual The residual function. + * @param numDofs The number of degrees of freedom. */ gsStaticOpt( const Residual_t &Residual, const index_t numDofs ) : @@ -131,9 +173,10 @@ class gsStaticOpt : public gsStaticBase } /** - * @brief Constructs a new instance. + * @brief Constructor for gsStaticOpt with arc-length residual function. * - * @param[in] ALResidual The residual as arc-length object + * @param ALResidual The arc-length residual function. + * @param numDofs The number of degrees of freedom. */ gsStaticOpt( const ALResidual_t & ALResidual, const index_t numDofs ) : @@ -146,6 +189,11 @@ class gsStaticOpt : public gsStaticBase } public: + /** + * @brief Returns the optimizer options. + * + * @return A reference to the optimizer options. + */ gsOptionList & optimizerOptions() { return m_optimizer.options(); } /// gsStaticBase base functions @@ -158,7 +206,7 @@ class gsStaticOpt : public gsStaticBase // /// See \ref gsStaticBase // void initOutput() override; - + // /// See \ref gsStaticBase // void stepOutput(index_t k) override; @@ -183,22 +231,12 @@ class gsStaticOpt : public gsStaticBase GISMO_NO_IMPLEMENTATION; } -// Class-specific functions -// public: -// void initialize(); - protected: /// See \ref solve() void _solve(); /// Initializes the method void _init(); -// gsVector _computeResidual(const gsVector & U); - -// public: - -// /// Return the residual norm -// T residualNorm() const { return m_R.norm(); } protected: Optimizer m_optimizer; diff --git a/src/gsStaticSolvers/gsStaticOpt.hpp b/src/gsStaticSolvers/gsStaticOpt.hpp index 7ca2c6a..101e0ea 100644 --- a/src/gsStaticSolvers/gsStaticOpt.hpp +++ b/src/gsStaticSolvers/gsStaticOpt.hpp @@ -93,150 +93,11 @@ gsStatus gsStaticOpt::solve() { this->getOptions(); m_optimizer.solve(m_U+m_DeltaU); - + m_DeltaU = m_optimizer.currentDesign() - m_U; m_U += m_DeltaU; m_numIterations = m_optimizer.iterations(); return gsStatus::Success; } - - - -//! [OptProblem] - -// template -// void gsStaticOpt::getOptions() -// { -// Base::getOptions(); -// } - -// template -// gsStatus gsStaticOpt::solve() -// { -// try -// { -// _solve(); -// 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 gsStaticOpt::_solve() -// { -// m_Eks.clear(); -// m_Eks.reserve(m_maxIterations); - -// if (m_verbose) initOutput(); - -// _start(); - -// m_Ek0 = m_Ek; -// m_Eks.push_back(m_Ek); - -// if (m_verbose != 0) stepOutput(0); -// index_t resetIt = 0; -// for (m_numIterations=1; m_numIterations!=m_maxIterations; m_numIterations++, resetIt++) -// { -// _iteration(); -// if ((m_c==0 && m_Ek_prev > m_Ek) || resetIt==m_resetIterations)// || (m_Ek/m_Ek_prev > 1/m_tolE && m_Ek_prev!=0)) -// { -// resetIt = 0; -// _peak(); -// } - -// if (m_verbose!=0) -// if (m_numIterations % m_verbose == 0 || m_verbose==-1 ) stepOutput(m_numIterations); - -// m_Eks.push_back(m_Ek); - -// m_residualOld = m_residual; - -// if (m_residual/m_residualIni < m_tolF && m_Ek/m_Ek0 < m_tolE && m_deltaU.norm()/m_DeltaU.norm() < m_tolU) -// { -// m_U += m_DeltaU; -// gsDebug <<"Converged: \n"; -// gsDebug <<"\t |R|/|R0| = "< -// void gsStaticOpt::initialize() -// { -// this->reset(); -// getOptions(); -// } - -// template -// void gsStaticOpt::reset() -// { -// Base::reset(); -// m_status = gsStatus::NotStarted; -// } - -// template -// void gsStaticOpt::_init() -// { -// this->reset(); -// defaultOptions(); -// } - -// template -// gsVector gsStaticOpt::_computeResidual(const gsVector & U) -// { -// gsVector resVec; -// if (!m_residualFun(U, resVec)) -// throw 2; -// return resVec; -// } - -// template -// gsOptProblemStatic::gsOptProblemStatic() -// { - -// // Number of design variables -// m_numDesignVars = 2; -// // Number of constraints -// m_numConstraints = 0; -// // Number of non-zeros in the Jacobian of the constraints -// m_numConJacNonZero = 0; - -// m_desLowerBounds.resize(m_numDesignVars); -// m_desUpperBounds.resize(m_numDesignVars); - -// m_desLowerBounds.setConstant(-1e19); -// m_desUpperBounds.setConstant( 1e19); - -// // we initialize x in bounds, in the upper right quadrant -// m_curDesign.resize(m_numDesignVars,1); -// m_curDesign.setZero(); -// } - } // namespace gismo diff --git a/tutorials/CMakeLists.txt b/tutorials/CMakeLists.txt index 54267b8..a5d6550 100644 --- a/tutorials/CMakeLists.txt +++ b/tutorials/CMakeLists.txt @@ -11,7 +11,7 @@ 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}") + set_property(TEST ${tarname} PROPERTY LABELS "${PROJECT_NAME}-tutorials") if(GISMO_BUILD_EXAMPLES) set_target_properties(${tarname} PROPERTIES FOLDER "${PROJECT_NAME}") else(GISMO_BUILD_EXAMPLES) diff --git a/unittests/gsStaticSolver_test.cpp b/unittests/gsStaticSolver_test.cpp new file mode 100644 index 0000000..6cb0378 --- /dev/null +++ b/unittests/gsStaticSolver_test.cpp @@ -0,0 +1,421 @@ +/** @file gsThinShellAssembler_test.cpp + + @brief Provides unittests for the gsThinShellAssembler class + + * Balloon: unit-test based on a hyperelastic balloon inflated by a follower pressure. + This test allows to test the follower pressure as well as the hyperelastic material models (incompressible) + + * UAT: unit-test based on a uni-axial tension test + This test allows to test the hyperelastic material models (incompressible and compressible) + + * Modal: unit-test based on a modal analysis + This test allows to test the mass matrix and the compressible material model + + + == BASIC REFERENCE == + - TEST(NAME_OF_TEST) { body_of_test } + - TEST_FIXTURE(NAME_OF_FIXTURE,NAME_OF_TEST){ body_of_test } + + == CHECK MACRO REFERENCE == + - CHECK(EXPR); + - CHECK_EQUAL(EXPECTED,ACTUAL); + - CHECK_CLOSE(EXPECTED,ACTUAL,EPSILON); + - CHECK_ARRAY_EQUAL(EXPECTED,ACTUAL,LENGTH); + - CHECK_ARRAY_CLOSE(EXPECTED,ACTUAL,LENGTH,EPSILON); + - CHECK_ARRAY2D_EQUAL(EXPECTED,ACTUAL,ROWCOUNT,COLCOUNT); + - CHECK_ARRAY2D_CLOSE(EXPECTED,ACTUAL,ROWCOUNT,COLCOUNT,EPSILON); + - CHECK_THROW(EXPR,EXCEPTION_TYPE_EXPECTED); + + == TIME CONSTRAINTS == + - UNITTEST_TIME_CONSTRAINT(TIME_IN_MILLISECONDS); + - UNITTEST_TIME_CONSTRAINT_EXEMPT(); + + == MORE INFO == + See: https://unittest-cpp.github.io/ + + Author(s): H.M.Verhelst (2019 - ..., TU Delft) + **/ + +#include "gismo_unittest.h" // Brings in G+Smo and the UnitTest++ framework +#ifdef gsKLShell_ENABLED +#include +#endif + +#ifdef gsOptim_ENABLED +#include +#endif + +#include +#include +#include +#include + +SUITE(gsThinShellAssembler_test) // The suite should have the same name as the file +{ + +#ifdef gsKLShell_ENABLED + // solver: 0: newton, 1: DR, 2: OPT + std::pair UAT_analytical(const std::vector solver, const index_t material=1, const index_t impl=1, const bool Compressibility=false); + std::pair UAT_numerical(const std::vector solver, const index_t material=1, const index_t impl=1, const bool Compressibility=false); + void UAT_CHECK(const std::vector solver, const index_t material=1, const index_t impl=1, const bool Compressibility=false); +#endif + +///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +#ifdef gsKLShell_ENABLED + TEST(StaticSolver_UAT_NR) + { + std::vector solver({0}); + UAT_CHECK(solver); + } + + // TEST(StaticSolver_UAT_DR) + // { + // std::vector solver({1}); + // UAT_CHECK(solver); + // } + +#ifdef gsOptim_ENABLED + TEST(StaticSolver_UAT_OP) + { + std::vector solver({2}); + UAT_CHECK(solver); + } +#endif + +#ifdef gsOptim_ENABLED + TEST(StaticSolver_UAT_OP_NR) + { + std::vector solver({2,0}); + UAT_CHECK(solver); + } +#endif + + TEST(StaticSolver_UAT_DR_NR) + { + std::vector solver({1,0}); + UAT_CHECK(solver); + } + +#ifdef gsOptim_ENABLED + TEST(StaticSolver_UAT_DR_OP) + { + std::vector solver({1,2}); + UAT_CHECK(solver); + } +#endif + + std::pair UAT_numerical(const std::vector solver, const index_t material, const index_t impl, const bool Compressibility) + { + //! [Parse command line] + index_t numRefine = 1; + index_t numElevate = 1; + + real_t E_modulus = 1.0; + real_t PoissonRatio; + real_t Density = 1.0; + real_t Ratio = 7.0; + + real_t mu = 1.5e6; + real_t thickness = 0.001; + + real_t alpha1,alpha2,alpha3,mu1,mu2,mu3; + alpha1 = 1.3; + mu1 = 6.3e5/4.225e5*mu; + alpha2 = 5.0; + mu2 = 0.012e5/4.225e5*mu; + alpha3 = -2.0; + mu3 = -0.1e5/4.225e5*mu; + + if (!Compressibility) + PoissonRatio = 0.5; + else + PoissonRatio = 0.45; + + E_modulus = 2*mu*(1+PoissonRatio); + + //! [Parse command line] + + //! [Read input file] + gsMultiPatch<> mp, mp_def; + + mp.addPatch( gsNurbsCreator<>::BSplineSquare(1) ); // degree + + if (numElevate!=0) + mp.degreeElevate(numElevate); + + // h-refine + for (int r =0; r < numRefine; ++r) + mp.uniformRefine(); + + mp_def = mp; + + //! [Refinement] + gsMultiBasis<> dbasis(mp); + + gsBoundaryConditions<> bc; + bc.setGeoMap(mp); + + gsPointLoads pLoads = gsPointLoads(); + + real_t lambda = 2.0; + gsConstantFunction<> displx(lambda-1.0,2); + + GISMO_ENSURE(mp.targetDim()==2,"Geometry must be planar (targetDim=2)!"); + bc.addCondition(boundary::west, condition_type::dirichlet, 0, 0, false, 0 ); + + bc.addCondition(boundary::east, condition_type::dirichlet, &displx, 0, false, 0 ); + + bc.addCondition(boundary::south, condition_type::dirichlet, 0, 0, false, 1 ); + + //! [Refinement] + + // Linear isotropic material model + gsVector<> tmp(2); + tmp.setZero(); + gsConstantFunction<> force(tmp,2); + gsFunctionExpr<> t(std::to_string(thickness),2); + gsFunctionExpr<> E(std::to_string(E_modulus),2); + gsFunctionExpr<> nu(std::to_string(PoissonRatio),2); + gsFunctionExpr<> rho(std::to_string(Density),2); + gsConstantFunction<> ratio(Ratio,2); + + gsConstantFunction<> alpha1fun(alpha1,2); + gsConstantFunction<> mu1fun(mu1,2); + gsConstantFunction<> alpha2fun(alpha2,2); + gsConstantFunction<> mu2fun(mu2,2); + gsConstantFunction<> alpha3fun(alpha3,2); + gsConstantFunction<> mu3fun(mu3,2); + + std::vector*> parameters(3); + parameters[0] = &E; + parameters[1] = ν + parameters[2] = ∶ + gsMaterialMatrixBase* materialMatrix; + + if (material==4) + { + parameters.resize(8); + parameters[0] = &E; + parameters[1] = ν + parameters[2] = &mu1fun; + parameters[3] = &alpha1fun; + parameters[4] = &mu2fun; + parameters[5] = &alpha2fun; + parameters[6] = &mu3fun; + parameters[7] = &alpha3fun; + } + + gsOptionList options; + if (material==0) + { + GISMO_ERROR("This test is not available for SvK models"); + } + else + { + options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",material); + options.addSwitch("Compressibility","Compressibility: (false): Imcompressible | (true): Compressible",Compressibility); + options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",impl); + materialMatrix = getMaterialMatrix<2,real_t>(mp,t,parameters,rho,options); + } + + gsThinShellAssemblerBase* assembler; + assembler = new gsThinShellAssembler<2, real_t, false >(mp,dbasis,bc,force,materialMatrix); + + assembler->setPointLoads(pLoads); + + // Function for the Jacobian + typedef std::function (gsVector const &)> Jacobian_t; + typedef std::function (gsVector const &) > Residual_t; + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&mp_def](gsVector const &x, gsSparseMatrix & m) + { + ThinShellAssemblerStatus status; + assembler->constructSolution(x,mp_def); + status = assembler->assembleMatrix(mp_def); + m = assembler->matrix(); + return status == ThinShellAssemblerStatus::Success; + }; + + gsStructuralAnalysisOps::Residual_t Residual = [&assembler,&mp_def](gsVector const &x, gsVector & result) + { + ThinShellAssemblerStatus status; + assembler->constructSolution(x,mp_def); + status = assembler->assembleVector(mp_def); + result = assembler->rhs(); + return status == ThinShellAssemblerStatus::Success; + }; + + // Define Matrices + assembler->assemble(); + gsSparseMatrix<> K = assembler->matrix(); + gsVector<> F = assembler->rhs(); + assembler->assembleMass(true); + gsVector<> M = assembler->rhs(); + + gsStaticDR DRM(M,F,Residual); + gsOptionList DROptions = DRM.options(); + DROptions.setReal("damping",0.0); + DROptions.setReal("alpha",1e9); + DROptions.setInt("maxIt",1e6); + DROptions.setReal("tolF",1e-8); + DROptions.setReal("tolU",1e-6); + DROptions.setReal("tolE",1e-4); + DROptions.setInt("verbose",0); + DRM.setOptions(DROptions); + DRM.initialize(); + +#ifdef gsOptim_ENABLED + gsStaticOpt::LBFGS> OPT(Residual,assembler->numDofs()); +#else + gsStaticOpt> OPT(Residual,assembler->numDofs()); +#endif + gsOptionList OPTOptions = OPT.options(); + OPTOptions.setInt("maxIt",100); + OPTOptions.setReal("tolF",1e-8); + OPTOptions.setReal("tolU",1e-6); + OPTOptions.setInt("verbose",0); + OPT.setOptions(OPTOptions); + OPT.initialize(); + + gsStaticNewton NWT(K,F,Jacobian,Residual); + gsOptionList NWTOptions = NWT.options(); + NWTOptions.setInt("maxIt",20); + NWTOptions.setReal("tolF",1e-8); + NWTOptions.setReal("tolU",1e-6); + NWTOptions.setInt("verbose",0); + NWT.setOptions(NWTOptions); + NWT.initialize(); + + std::vector*> solverVec(solver.size()); + for (size_t s=0; s!=solver.size(); s++) + if (solver[s]==0) + solverVec[s] = &NWT; + else if (solver[s]==1) + solverVec[s] = &DRM; + else if (solver[s]==2) + solverVec[s] = &OPT; + + gsStaticComposite compositeSolver(solverVec); + compositeSolver.initialize(); + compositeSolver.solve(); + CHECK(compositeSolver.status() == gsStatus::Success); + mp_def = assembler->constructSolution(compositeSolver.solution()); + + gsMultiPatch<> deformation = mp_def; + for (size_t k = 0; k != mp_def.nPatches(); ++k) + deformation.patch(k).coefs() -= mp.patch(k).coefs(); + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // Check solutions + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + // NOTE: all analytical solutions for compressible materials are fixed for displ=1; (lambda=2) + + // Compute stretches (should be the same everywhere) + // Ordering: lambda(0) < lambda(1); lambda(2) is ALWAYS the through-thickness stretch + gsVector<> pt(2); + pt<<1,0; + gsMatrix<> lambdas = assembler->computePrincipalStretches(pt,mp_def,0); + + // Get the total force on the tension boundary + patchSide ps(0,boundary::east); + gsMatrix<> forceVector = assembler->boundaryForce(mp_def,ps); + real_t sideForce = forceVector.sum(); + real_t S = -sideForce / (thickness*lambdas(0)*lambdas(2)); + real_t L = lambdas(0); + + std::pair result; + result.first = L; + result.second = S; + + delete materialMatrix; + delete assembler; + + return result; + } + + std::pair UAT_analytical(const std::vector solver, const index_t material, const index_t impl, const bool Compressibility) + { + real_t PoissonRatio; + real_t Ratio = 7.0; + + real_t mu = 1.5e6; + + real_t alpha1,alpha2,alpha3,mu1,mu2,mu3; + alpha1 = 1.3; + mu1 = 6.3e5/4.225e5*mu; + alpha2 = 5.0; + mu2 = 0.012e5/4.225e5*mu; + alpha3 = -2.0; + mu3 = -0.1e5/4.225e5*mu; + + if (!Compressibility) + PoissonRatio = 0.5; + else + PoissonRatio = 0.45; + + real_t lambda = 2.0; + + real_t San,J,K,Lan; + if (material==1 && Compressibility) + { + K = 2*mu*(1+PoissonRatio)/(3-6*PoissonRatio); + J = 1.105598565;// specific for lambda==2!! + San = lambda*(0.5*mu*(-(2*(math::pow(lambda,2)+2*J/lambda))/(3*math::pow(J,2./3.)*lambda)+2*lambda/math::pow(J,2./3.))+0.25*K*(2*math::pow(J,2)/lambda-2./lambda))/J; + Lan = math::pow(J/lambda,0.5); + } + else if (material==1 && !Compressibility) + { + San = mu * (lambda*lambda - 1/lambda); + Lan = math::pow(1./lambda,0.5); + } + else if (material==3 && Compressibility) + { + real_t c2 = 1.0 / (Ratio+1); + real_t c1 = 1.0 - c2; + K = 2*mu*(1+PoissonRatio)/(3-6*PoissonRatio); + J = 1.099905842;// specific for lambda==2!! + San = lambda*(0.5*c1*mu*(-(2*(math::pow(lambda,2)+2*J/lambda))/(3*math::pow(J,2./3.)*lambda)+2*lambda/math::pow(J,2./3.))+0.5*c2*mu*(-(4*(2*lambda*J+math::pow(J,2)/math::pow(lambda,2)))/(3*math::pow(J,4./3.)*lambda)+4/math::pow(J,1./3.))+0.25*K*(2*math::pow(J,2)/lambda-2/lambda))/J; + Lan = math::pow(J/lambda,0.5); + } + else if (material==3 && !Compressibility) + { + real_t c2 = 1.0 / (Ratio+1); + real_t c1 = 1.0 - c2; + San =-mu*(c2*lambda*lambda+c2/lambda+c1)/lambda+lambda*(c1*lambda*mu+2*c2*mu); + Lan = math::pow(1./lambda,0.5); + } + else if (material==4 && Compressibility) + { + K = 2*mu*(1+PoissonRatio)/(3-6*PoissonRatio); + J = 1.088778638;// specific for lambda==2!! + San = 1./J* (lambda *( mu1*(2*math::pow(lambda/math::pow(J,1./3.),alpha1)*alpha1/(3*lambda)-2*math::pow(math::pow(J/lambda,0.5)/math::pow(J,1./3.),alpha1)*alpha1/(3*lambda))/alpha1+mu2*(2*math::pow(lambda/math::pow(J,1./3.),alpha2)*alpha2/(3*lambda)-2*math::pow(math::pow(J/lambda,0.5)/math::pow(J,1./3.),alpha2)*alpha2/(3*lambda))/alpha2+mu3*(2*math::pow(lambda/math::pow(J,1./3.),alpha3)*alpha3/(3*lambda)-2*math::pow(math::pow(J/lambda,0.5)/math::pow(J,1./3.),alpha3)*alpha3/(3*lambda))/alpha3+0.25*K*(2*math::pow(J,2)/lambda-2/lambda) ) ); + Lan = math::pow(J/lambda,0.5); + } + else if (material==4 && !Compressibility) + { + San =-mu1*math::pow((1./lambda),0.5*alpha1)-mu2*math::pow((1./lambda),0.5*alpha2)-mu3*math::pow((1./lambda),0.5*alpha3)+mu1*math::pow(lambda,alpha1)+mu2*math::pow(lambda,alpha2)+mu3*math::pow(lambda,alpha3); + Lan = math::pow(1./lambda,0.5); + } + else + GISMO_ERROR("Material not treated"); + + std::pair result; + result.first = Lan; + result.second = San; + return result; + } + + void UAT_CHECK(const std::vector solver, const index_t material, const index_t impl, const bool Compressibility) + { + if (material==4 && impl!=3) + CHECK(true); + + real_t Lnum, Snum, Lana, Sana; + std::tie(Lnum,Snum) = UAT_numerical(solver,material,impl,Compressibility); + std::tie(Lana,Sana) = UAT_analytical(solver,material,impl,Compressibility); + CHECK_CLOSE(std::abs(Lnum-Lana)/Lana,0,1e-7); + } +#endif + +} From 25b6263aa55c44966e1fc68b5fc65b5a41ec8226 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Sun, 8 Dec 2024 19:51:10 +0100 Subject: [PATCH 7/9] change optim to HLBFGS --- .github/workflows/ci.yml | 2 +- unittests/gsStaticSolver_test.cpp | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 772260b..429d352 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -49,4 +49,4 @@ jobs: - name: "Run for ${{ matrix.os }}" shell: bash working-directory: ${{runner.workspace}} - run: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" -D CTEST_CONFIGURATION_TYPE=RelWithDebInfo -D LABELS_FOR_SUBPROJECTS="gsStructuralAnalysis-benchmarks" -D CTEST_SITE="${{ matrix.os }}_[actions]" -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON" -D GISMO_OPTIONAL='gsOptim\\;gsSpectra\\;gsKLShell\\;gsElasticity\\;${{ github.event.repository.name }}' -Q \ No newline at end of file + run: ctest -S ${{ github.event.repository.name }}/gismo/cmake/ctest_script.cmake -D CTEST_BUILD_NAME="${{ github.event.repository.name }}_actions_$GITHUB_RUN_NUMBER" -D CTEST_CONFIGURATION_TYPE=RelWithDebInfo -D LABELS_FOR_SUBPROJECTS="gsStructuralAnalysis-benchmarks" -D CTEST_SITE="${{ matrix.os }}_[actions]" -D CMAKE_ARGS="-DCMAKE_BUILD_TYPE=$BUILD_TYPE;-DCMAKE_CXX_STANDARD=11;-DGISMO_WITH_XDEBUG=ON;-DGISMO_BUILD_UNITTESTS=ON" -D GISMO_OPTIONAL='gsHLBFGS\\;gsSpectra\\;gsKLShell\\;gsElasticity\\;${{ github.event.repository.name }}' -Q \ No newline at end of file diff --git a/unittests/gsStaticSolver_test.cpp b/unittests/gsStaticSolver_test.cpp index 6cb0378..33e936a 100644 --- a/unittests/gsStaticSolver_test.cpp +++ b/unittests/gsStaticSolver_test.cpp @@ -41,8 +41,8 @@ #include #endif -#ifdef gsOptim_ENABLED -#include +#ifdef gsHLBFGS_ENABLED +#include #endif #include @@ -75,7 +75,7 @@ SUITE(gsThinShellAssembler_test) // The suite should have the sa // UAT_CHECK(solver); // } -#ifdef gsOptim_ENABLED +#ifdef gsHLBFGS_ENABLED TEST(StaticSolver_UAT_OP) { std::vector solver({2}); @@ -83,7 +83,7 @@ SUITE(gsThinShellAssembler_test) // The suite should have the sa } #endif -#ifdef gsOptim_ENABLED +#ifdef gsHLBFGS_ENABLED TEST(StaticSolver_UAT_OP_NR) { std::vector solver({2,0}); @@ -97,7 +97,7 @@ SUITE(gsThinShellAssembler_test) // The suite should have the sa UAT_CHECK(solver); } -#ifdef gsOptim_ENABLED +#ifdef gsHLBFGS_ENABLED TEST(StaticSolver_UAT_DR_OP) { std::vector solver({1,2}); @@ -265,8 +265,8 @@ SUITE(gsThinShellAssembler_test) // The suite should have the sa DRM.setOptions(DROptions); DRM.initialize(); -#ifdef gsOptim_ENABLED - gsStaticOpt::LBFGS> OPT(Residual,assembler->numDofs()); +#ifdef gsHLBFGS_ENABLED + gsStaticOpt> OPT(Residual,assembler->numDofs()); #else gsStaticOpt> OPT(Residual,assembler->numDofs()); #endif From 2fcf61421a9398fa5fc42aea52daa75d2d480081 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Tue, 10 Dec 2024 23:16:29 +0100 Subject: [PATCH 8/9] fix unittest CI --- unittests/gsStaticSolver_test.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/unittests/gsStaticSolver_test.cpp b/unittests/gsStaticSolver_test.cpp index 33e936a..95328f8 100644 --- a/unittests/gsStaticSolver_test.cpp +++ b/unittests/gsStaticSolver_test.cpp @@ -191,7 +191,7 @@ SUITE(gsThinShellAssembler_test) // The suite should have the sa parameters[0] = &E; parameters[1] = ν parameters[2] = ∶ - gsMaterialMatrixBase* materialMatrix; + gsMaterialMatrixBase::uPtr materialMatrix; if (material==4) { @@ -328,7 +328,6 @@ SUITE(gsThinShellAssembler_test) // The suite should have the sa result.first = L; result.second = S; - delete materialMatrix; delete assembler; return result; From 5cffa404ad43c6e54cd22e00556e2bd7785b4d5f Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Tue, 10 Dec 2024 23:38:58 +0100 Subject: [PATCH 9/9] unittest bug --- unittests/gsStaticSolver_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unittests/gsStaticSolver_test.cpp b/unittests/gsStaticSolver_test.cpp index 95328f8..bff97e0 100644 --- a/unittests/gsStaticSolver_test.cpp +++ b/unittests/gsStaticSolver_test.cpp @@ -50,7 +50,7 @@ #include #include -SUITE(gsThinShellAssembler_test) // The suite should have the same name as the file +SUITE(gsStaticSolver_test) // The suite should have the same name as the file { #ifdef gsKLShell_ENABLED