diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 25efed2..d7e99cf 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -37,4 +37,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="gsElasticity-examples" -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="${{ 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="gsElasticity-examples;gsElasticity-tutorials" -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="${{ github.event.repository.name }}" -Q diff --git a/CMakeLists.txt b/CMakeLists.txt index 8cfdd48..392354d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,11 +55,15 @@ add_definitions(-DELAST_DATA_DIR="${CMAKE_CURRENT_SOURCE_DIR}/filedata/") # add example files if(GISMO_BUILD_EXAMPLES) add_custom_target(${PROJECT_NAME}-examples) + add_custom_target(${PROJECT_NAME}-tutorials) add_subdirectory(examples) + add_subdirectory(tutorials) get_property(dirs TARGET ${PROJECT_NAME}-examples PROPERTY BUILDSYSTEM_TARGETS) + get_property(dirs TARGET ${PROJECT_NAME}-tutorials PROPERTY BUILDSYSTEM_TARGETS) message(STATUS "dirs: ${dirs}") else() add_subdirectory(examples EXCLUDE_FROM_ALL) + add_subdirectory(tutorials EXCLUDE_FROM_ALL) endif(GISMO_BUILD_EXAMPLES) diff --git a/filedata/paraboloid_volume.xml b/filedata/paraboloid_volume.xml new file mode 100644 index 0000000..38da6e8 --- /dev/null +++ b/filedata/paraboloid_volume.xml @@ -0,0 +1,55 @@ + + + + + + + 0 0 0 1 1 1 + + + 0 0 0 1 1 1 + + + 0 0 0 1 1 1 + + + + 0 0 0.49 + 0.5 0 0.99 + 1 0 0.49 + 0 0.5 0.99 + 0.5 0.5 1.45 + 1 0.5 0.99 + 0 1 0.49 + 0.5 1 0.99 + 1 1 0.49 + 0 0 0.5 + 0.5 0 1 + 1 0 0.5 + 0 0.5 1 + 0.5 0.5 1.5 + 1 0.5 1 + 0 1 0.5 + 0.5 1 1 + 1 1 0.5 + 0 0 0.51 + 0.5 0 1.01 + 1 0 0.51 + 0 0.5 1.01 + 0.5 0.5 1.51 + 1 0.5 1.01 + 0 1 0.51 + 0.5 1 1.01 + 1 1 0.51 + + + + 0 0 + 0 1 +0 2 +0 3 +0 4 + + + + diff --git a/gsElasticityAssembler.hpp b/gsElasticityAssembler.hpp index 78c867e..6e24899 100644 --- a/gsElasticityAssembler.hpp +++ b/gsElasticityAssembler.hpp @@ -146,7 +146,7 @@ void gsElasticityAssembler::assemble(bool saveEliminationMatrix) // Compute volumetric integrals and write to the global linear system if (m_bases.size() == unsigned(m_dim)) // displacement formulation { - GISMO_ENSURE(m_options.getInt("MaterialLaw") == material_law::hooke, + GISMO_ENSURE(m_options.getInt("MaterialLaw") == material_law::hooke || m_options.getInt("MaterialLaw") == material_law::saint_venant_kirchhoff, "Material law not specified OR not supported!"); if (saveEliminationMatrix) { @@ -189,7 +189,7 @@ bool gsElasticityAssembler::assemble(const gsMatrix & solutionVector, if (checkDisplacement(m_pde_ptr->patches(),displacement) != -1) return false; - if (m_bases.size() == unsigned(m_dim)) // displacement formulation + if (m_bases.size() == unsigned(m_dim)) // displacement formulation assemble(displacement); else // mixed formulation (displacement + pressure) { diff --git a/tutorials/CMakeLists.txt b/tutorials/CMakeLists.txt new file mode 100644 index 0000000..18a0414 --- /dev/null +++ b/tutorials/CMakeLists.txt @@ -0,0 +1,25 @@ +###################################################################### +## CMakeLists.txt --- gsStructuralAnalysis/tutorials +## 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}") + add_dependencies(${PROJECT_NAME}-tutorials ${tarname}) + else(GISMO_BUILD_EXAMPLES) + set_target_properties(${tarname} PROPERTIES + FOLDER "${PROJECT_NAME}" + EXCLUDE_FROM_ALL TRUE) + endif(GISMO_BUILD_EXAMPLES) + # 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/linear_solid.cpp b/tutorials/linear_solid.cpp new file mode 100644 index 0000000..c5c3bf9 --- /dev/null +++ b/tutorials/linear_solid.cpp @@ -0,0 +1,153 @@ +/** @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] +#include +//! [Includes] + +using namespace gismo; + +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","-1e4",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("DirichletValues",dirichlet::l2Projection); + + //! [Assemble linear part] + gsInfo<<"Assembling...\n"; + gsStopwatch clock; + clock.restart(); + assembler.assemble(); + gsInfo << "Assembled a system with " + << assembler.numDofs() << " dofs in " << clock.stop() << "s.\n"; + + gsInfo << "Solving...\n"; + clock.restart(); + + gsSparseMatrix<> matrix = assembler.matrix(); + gsVector<> vector = assembler.rhs(); + //! [Assemble linear part] + + //! [Solve linear problem] + gsInfo<<"Solving system with "< solVector; + gsSparseSolver<>::CGDiagonal solver; + solver.compute( matrix ); + solVector = solver.solve(vector); + //! [Solve linear 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 + // constructing stresses + 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; + +}// end main \ No newline at end of file