From 5b82d321f8b1c31b63c94c027694df7e657e26a2 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 8 Jan 2025 12:20:38 +0100 Subject: [PATCH] add tutorials for elasticity --- src/gsStaticSolvers/gsStaticBase.h | 7 +- tutorials/nonlinear_shell_static.cpp | 2 +- tutorials/nonlinear_solid_arcLength.cpp | 172 ++++++++++++++++++++ tutorials/nonlinear_solid_dynamic.cpp | 199 ++++++++++++++++++++++++ tutorials/nonlinear_solid_static.cpp | 181 +++++++++++++++++++++ 5 files changed, 557 insertions(+), 4 deletions(-) create mode 100644 tutorials/nonlinear_solid_arcLength.cpp create mode 100644 tutorials/nonlinear_solid_dynamic.cpp create mode 100644 tutorials/nonlinear_solid_static.cpp diff --git a/src/gsStaticSolvers/gsStaticBase.h b/src/gsStaticSolvers/gsStaticBase.h index 488b67b..85fefba 100644 --- a/src/gsStaticSolvers/gsStaticBase.h +++ b/src/gsStaticSolvers/gsStaticBase.h @@ -37,7 +37,7 @@ template class gsStaticBase { protected: - + typedef typename gsStructuralAnalysisOps::Residual_t Residual_t; typedef typename gsStructuralAnalysisOps::ALResidual_t ALResidual_t; typedef typename gsStructuralAnalysisOps::Jacobian_t Jacobian_t; @@ -89,7 +89,8 @@ class gsStaticBase virtual void setOptions(gsOptionList & options) {m_options.update(options,gsOptionList::addIfUnknown); } /// Get options - virtual gsOptionList options() const {return m_options;} + virtual const gsOptionList & options() const {return m_options;} + virtual gsOptionList & options() {return m_options;} /// Access the solution virtual gsVector solution() const {return m_U;} @@ -199,7 +200,7 @@ class gsStaticBase if (es.info()!=Spectra::CompInfo::Successful) { gsWarn<<"Spectra did not converge!\n"; - return false; + return false; } // if (es.info()==Spectra::CompInfo::NotComputed) diff --git a/tutorials/nonlinear_shell_static.cpp b/tutorials/nonlinear_shell_static.cpp index b4772b5..ab6111c 100644 --- a/tutorials/nonlinear_shell_static.cpp +++ b/tutorials/nonlinear_shell_static.cpp @@ -145,7 +145,7 @@ int main(int argc, char *argv[]) //! [Set static solver] gsStaticNewton solver(K,F,Jacobian,Residual); - solver.options().setInt("verbose",2); + solver.options().setInt("verbose",1); solver.initialize(); //! [Set static solver] diff --git a/tutorials/nonlinear_solid_arcLength.cpp b/tutorials/nonlinear_solid_arcLength.cpp new file mode 100644 index 0000000..0119373 --- /dev/null +++ b/tutorials/nonlinear_solid_arcLength.cpp @@ -0,0 +1,172 @@ +/** @file gsKLShell/tutorials/linear_solid.cpp + + @brief Tutorial for assembling solid elasticity + + 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 gsElasticity_ENABLED +#include +#include +#include +#endif + + +#include +#include +//! [Includes] + +using namespace gismo; + +#ifdef gsElasticity_ENABLED +int main(int argc, char *argv[]) +{ + //! [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; + std::string fileName = "paraboloid_volume.xml"; + gsReadFile<>(fileName, ori); + //! [Read geometry] + + //! [Initialize geometry] + // p-refine + if (numElevate!=0) + ori.degreeElevate(numElevate); + // h-refine (only in first two directions) + for (int r =0; r < numRefine; ++r) + { + ori.uniformRefine(1,1,0); + ori.uniformRefine(1,1,1); + } + + // creating basis + gsMultiBasis<> basis(ori); + //! [Initialize geometry] + + //! [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 + gsFunctionExpr<> surf_force("0","0","-1e10",3); + for (index_t c=0; c!=3; c++) + { + bc.addCornerValue(boundary::southwestfront, 0.0, 0, c); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::southeastfront, 0.0, 0, c); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northwestfront, 0.0, 0, c); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northeastfront, 0.0, 0, c); // (corner,value, patch, unknown) + } + bc.addCondition(boundary::front, condition_type::neumann, &surf_force); + //! [Set boundary conditions] + + //! [Set surface force] + // The surface force is defined in the physical space, i.e. 3D + gsFunctionExpr<> body_force("0","0","0",3); + //! [Set surface force] + + //! [Define the assembler] + gsElasticityAssembler assembler(ori,basis,bc,body_force); + assembler.options().setReal("YoungsModulus",1.E9); + assembler.options().setReal("PoissonsRatio",0.45); + assembler.options().setInt("MaterialLaw",material_law::saint_venant_kirchhoff); + assembler.options().setInt("DirichletValues",dirichlet::l2Projection); + + //! [Assemble linear part] + assembler.assemble(); + gsSparseMatrix<> K = assembler.matrix(); + gsVector<> F = assembler.rhs(); + //! [Assemble linear part] + + //! [Define nonlinear residual functions] + std::vector > fixedDofs = assembler.allFixedDofs(); + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&fixedDofs](gsVector const &x, gsSparseMatrix & m) + { + assembler.assemble(x,fixedDofs); + m = assembler.matrix(); + return true; + }; + + // Function for the Residual + gsStructuralAnalysisOps::ALResidual_t ALResidual = [&fixedDofs,&assembler,&F](gsVector const &x, real_t lam, gsVector & v) + { + assembler.assemble(x,fixedDofs); + v = F - lam * F - assembler.rhs(); // assembler rhs - force = Finternal + return true; + }; + //! [Define nonlinear residual functions] + + //! [Set ALM solver] + gsInfo<<"Solving system with "< solver(Jacobian,ALResidual,F); + solver.options().setSwitch("Verbose",true); + solver.setLength(0.1); + 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 displ; + assembler.constructSolution(solVector,assembler.allFixedDofs(),displ); + gsMultiPatch<> def = ori; + for (index_t p = 0; p < def.nPieces(); ++p) + def.patch(p).coefs() += displ.patch(p).coefs(); + + // ! [Export visualization in ParaView] + if (plot) + { + // Plot the displacements on the deformed geometry + gsField<> solField(def, displ); + std::string outputName = "Deformation" + std::to_string(k) + "_"; + gsWriteParaview<>( solField, outputName, 1000, true); + collection.addPart(outputName + "0.vts",k); + } + } + //! [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 gsElasticity enabled"); + return EXIT_FAILED; +#endif + +}// end main \ No newline at end of file diff --git a/tutorials/nonlinear_solid_dynamic.cpp b/tutorials/nonlinear_solid_dynamic.cpp new file mode 100644 index 0000000..cd9378a --- /dev/null +++ b/tutorials/nonlinear_solid_dynamic.cpp @@ -0,0 +1,199 @@ +/** @file gsKLShell/tutorials/linear_solid.cpp + + @brief Tutorial for assembling solid elasticity + + 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 gsElasticity_ENABLED +#include +#include +#include +#endif + + +#include +#include +//! [Includes] + +using namespace gismo; + +#ifdef gsElasticity_ENABLED +int main(int argc, char *argv[]) +{ + //! [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; + std::string fileName = "paraboloid_volume.xml"; + gsReadFile<>(fileName, ori); + //! [Read geometry] + + //! [Initialize geometry] + // p-refine + if (numElevate!=0) + ori.degreeElevate(numElevate); + // h-refine (only in first two directions) + for (int r =0; r < numRefine; ++r) + { + ori.uniformRefine(1,1,0); + ori.uniformRefine(1,1,1); + } + + // creating basis + gsMultiBasis<> basis(ori); + //! [Initialize geometry] + + //! [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 + gsFunctionExpr<> surf_force("0","0","-1e5",3); +for (index_t c=0; c!=3; c++) + { + bc.addCornerValue(boundary::southwestfront, 0.0, 0, c); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::southeastfront, 0.0, 0, c); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northwestfront, 0.0, 0, c); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northeastfront, 0.0, 0, c); // (corner,value, patch, unknown) + } + bc.addCondition(boundary::front, condition_type::neumann, &surf_force); + //! [Set boundary conditions] + + //! [Set surface force] + // The surface force is defined in the physical space, i.e. 3D + gsFunctionExpr<> body_force("0","0","0",3); + //! [Set surface force] + + //! [Define the assembler] + gsElasticityAssembler assembler(ori,basis,bc,body_force); + assembler.options().setReal("YoungsModulus",1.E9); + assembler.options().setReal("PoissonsRatio",0.45); + assembler.options().setInt("MaterialLaw",material_law::saint_venant_kirchhoff); + assembler.options().setInt("DirichletValues",dirichlet::l2Projection); + + gsMassAssembler assemblerMass(ori,basis,bc,body_force); + assemblerMass.options() = assembler.options(); + assemblerMass.options().addReal("Density","Density of the material",1.E8); + + //! [Assemble linear part] + assembler.assemble(); + gsSparseMatrix<> K = assembler.matrix(); + gsVector<> F = assembler.rhs(); + + assemblerMass.assemble(); + gsSparseMatrix<> M = assemblerMass.matrix(); + //! [Assemble linear part] + + //! [Define nonlinear residual functions] + std::vector > fixedDofs = assembler.allFixedDofs(); + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&fixedDofs](gsVector const &x, gsSparseMatrix & m) + { + assembler.assemble(x,fixedDofs); + m = assembler.matrix(); + return true; + }; + + // Function for the Residual + gsStructuralAnalysisOps::Residual_t Residual = [&fixedDofs,&assembler](gsVector const &x, gsVector & v) + { + assembler.assemble(x,fixedDofs); + v = assembler.rhs(); + return true; + }; + //! [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 displ; + assembler.constructSolution(U,assembler.allFixedDofs(),displ); + gsMultiPatch<> def = ori; + for (index_t p = 0; p < def.nPieces(); ++p) + def.patch(p).coefs() += displ.patch(p).coefs(); + + // ! [Export visualization in ParaView] + if (plot) + { + // Plot the displacements on the deformed geometry + gsField<> 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 gsElasticity enabled"); + return EXIT_FAILED; +#endif + +}// end main \ No newline at end of file diff --git a/tutorials/nonlinear_solid_static.cpp b/tutorials/nonlinear_solid_static.cpp new file mode 100644 index 0000000..7529631 --- /dev/null +++ b/tutorials/nonlinear_solid_static.cpp @@ -0,0 +1,181 @@ +/** @file gsKLShell/tutorials/linear_solid.cpp + + @brief Tutorial for assembling solid elasticity + + 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 gsElasticity_ENABLED +#include +#include +#include +#endif + + +#include +#include +//! [Includes] + +using namespace gismo; + +#ifdef gsElasticity_ENABLED +int main(int argc, char *argv[]) +{ + //! [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; + std::string fileName = "paraboloid_volume.xml"; + gsReadFile<>(fileName, ori); + //! [Read geometry] + + //! [Initialize geometry] + // p-refine + if (numElevate!=0) + ori.degreeElevate(numElevate); + // h-refine (only in first two directions) + for (int r =0; r < numRefine; ++r) + { + ori.uniformRefine(1,1,0); + ori.uniformRefine(1,1,1); + } + + // creating basis + gsMultiBasis<> basis(ori); + //! [Initialize geometry] + + //! [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 + gsFunctionExpr<> surf_force("0","0","-1e5",3); + for (index_t c=0; c!=3; c++) + { + bc.addCornerValue(boundary::southwestfront, 0.0, 0, c); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::southeastfront, 0.0, 0, c); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northwestfront, 0.0, 0, c); // (corner,value, patch, unknown) + bc.addCornerValue(boundary::northeastfront, 0.0, 0, c); // (corner,value, patch, unknown) + } + bc.addCondition(boundary::front, condition_type::neumann, &surf_force); + bc.addCondition(boundary::front, condition_type::neumann, &surf_force); + //! [Set boundary conditions] + + //! [Set surface force] + // The surface force is defined in the physical space, i.e. 3D + gsFunctionExpr<> body_force("0","0","0",3); + //! [Set surface force] + + //! [Define the assembler] + gsElasticityAssembler assembler(ori,basis,bc,body_force); + assembler.options().setReal("YoungsModulus",1.E9); + assembler.options().setReal("PoissonsRatio",0.45); + assembler.options().setInt("MaterialLaw",material_law::saint_venant_kirchhoff); + assembler.options().setInt("DirichletValues",dirichlet::l2Projection); + + //! [Define nonlinear residual functions] + std::vector > fixedDofs = assembler.allFixedDofs(); + // Function for the Jacobian + gsStructuralAnalysisOps::Jacobian_t Jacobian = [&assembler,&fixedDofs](gsVector const &x, gsSparseMatrix & m) + { + assembler.assemble(x,fixedDofs); + m = assembler.matrix(); + return true; + }; + + // Function for the Residual + gsStructuralAnalysisOps::Residual_t Residual = [&fixedDofs,&assembler](gsVector const &x, gsVector & v) + { + assembler.assemble(x,fixedDofs); + v = assembler.rhs(); + return true; + }; + //! [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",1); + solver.initialize(); + //! [Set static solver] + + //! [Solve nonlinear problem] + gsInfo<<"Solving system with "< solVector = solver.solution(); + //! [Solve nonlinear problem] + + //! [Construct solution and deformed geometry] + // constructing solution as an IGA function + gsMultiPatch<> displ; + assembler.constructSolution(solVector,assembler.allFixedDofs(),displ); + gsMultiPatch<> def = ori; + for (index_t p = 0; p < def.nPieces(); ++p) + def.patch(p).coefs() += displ.patch(p).coefs(); + //! [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(3,1); + refPoint<<0.5,0.5,1.0; + // 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", 10000, true); + + // Plot the membrane Von-Mises stresses on the geometry + gsPiecewiseFunction<> VMm; + assembler.constructCauchyStresses(displ,VMm,stress_components::von_mises); + gsField<> membraneStress(def,VMm, true); + gsWriteParaview(membraneStress,"MembraneStress",10000); + } + // ! [Export visualization in ParaView] + + return EXIT_SUCCESS; +#else + GISMO_ERROR("The tutorial needs to be compiled with gsElasticity enabled"); + return EXIT_FAILED; +#endif + +}// end main \ No newline at end of file