diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..7465f11 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,114 @@ +name: CI + +on: + workflow_dispatch: + pull_request: + push: + +jobs: + build-and-test: + name: Run on ${{ matrix.os }} with SOFA ${{ matrix.sofa_branch }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + os: [ubuntu-22.04, macos-13, windows-2022] + sofa_branch: [master] + + steps: + + - name: Setup SOFA and environment + id: sofa + uses: sofa-framework/sofa-setup-action@v5 + with: + sofa_root: ${{ github.workspace }}/sofa + sofa_version: ${{ matrix.sofa_branch }} + sofa_scope: 'standard' + sofa_with_sofapython3: 'true' + + - name: Checkout source code + uses: actions/checkout@v2 + with: + path: ${{ env.WORKSPACE_SRC_PATH }} + + - name: Build and install + shell: bash + run: | + cmake_options="-GNinja \ + -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_INSTALL_PREFIX="$WORKSPACE_INSTALL_PATH" \ + -DCMAKE_PREFIX_PATH="$SOFA_ROOT/lib/cmake" \ + " + if [ -e "$(command -v ccache)" ]; then + cmake_options="$cmake_options -DCMAKE_C_COMPILER_LAUNCHER=ccache -DCMAKE_CXX_COMPILER_LAUNCHER=ccache" + fi + cmake_options="$(echo $cmake_options)" # prettify + + if [[ "$RUNNER_OS" == "Windows" ]]; then + cmd //c "${{ steps.sofa.outputs.vs_vsdevcmd }} \ + && cd /d $WORKSPACE_BUILD_PATH \ + && cmake $cmake_options ../src \ + && ninja install" + else + cd "$WORKSPACE_BUILD_PATH" + ccache -z + cmake $cmake_options ../src + ninja install + echo ${CCACHE_BASEDIR} + ccache -s + fi + + - name: Sanitize artifact name + id: sanitize + # This step removes special characters from the artifact name to ensure compatibility with upload-artifact + # Characters removed: " : < > | * ? \r \n \ / + # Spaces are replaced with underscores + # This sanitization prevents errors in artifact creation and retrieval + shell: pwsh + run: | + $originalName = "Shell_${{ steps.sofa.outputs.run_branch }}_for-SOFA-${{ steps.sofa.outputs.sofa_version }}_${{ runner.os }}" + $artifact_name = $originalName -replace '[":;<>|*?\r\n\\/]', '' -replace ' ', '_' + echo "artifact_name=$artifact_name" >> $env:GITHUB_OUTPUT + + - name: Create artifact + uses: actions/upload-artifact@v4.4.0 + with: + name: ${{ steps.sanitize.outputs.artifact_name }} + path: ${{ env.WORKSPACE_INSTALL_PATH }} + + - name: Install artifact + uses: actions/download-artifact@v4.1.7 + with: + name: ${{ steps.sanitize.outputs.artifact_name }} + path: ${{ env.WORKSPACE_ARTIFACT_PATH }} + + deploy: + name: Deploy artifacts + if: always() && startsWith(github.repository, 'SofaDefrost') && (startsWith(github.ref, 'refs/heads/') || startsWith(github.ref, 'refs/tags/')) # we are not on a fork and on a branch or a tag (not a PR) + needs: [build-and-test] + runs-on: ubuntu-latest + continue-on-error: true + steps: + - name: Get artifacts + uses: actions/download-artifact@v4.1.7 + with: + path: artifacts + + - name: Zip artifacts + shell: bash + run: | + cd $GITHUB_WORKSPACE/artifacts + for artifact in *; do + zip $artifact.zip -r $artifact/* + done + - name: Upload release + uses: softprops/action-gh-release@v1 + with: + name: ${{ github.ref_name }} + tag_name: release-${{ github.ref_name }} + fail_on_unmatched_files: false + target_commitish: ${{ github.ref_name }} + files: | + artifacts/Shell_*_Linux.zip + artifacts/Shell_*_Windows.zip + artifacts/Shell_*_macOS.zip diff --git a/CMakeLists.txt b/CMakeLists.txt index 5bdce1f..f4def9a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ cmake_minimum_required(VERSION 3.12) -project(SofaShells VERSION 1.0 LANGUAGES CXX) - +project(Shell VERSION 1.0 LANGUAGES CXX) +# Find and load CMake configuration of packages containing this plugin's dependencies find_package(Sofa.Config REQUIRED) sofa_find_package(Sofa.GL REQUIRED) sofa_find_package(Sofa.Type REQUIRED) @@ -17,113 +17,98 @@ sofa_find_package(Sofa.Component.Constraint.Lagrangian.Solver REQUIRED) sofa_find_package(Sofa.Component.Topology.Container.Dynamic REQUIRED) sofa_find_package(Sofa.Component.Collision.Detection.Intersection REQUIRED) - - -set(README_FILE "SofaShells.txt" ) +set(README_FILE README.md) set(SOFASHELLS_SRC_DIR src/SofaShells) - option(SOFA-PLUGIN_SHELLS_ADAPTIVITY "Enables shells adaptivity" OFF) - +# List all files +set(SHELL_SRC_DIR src/Shell) set(HEADER_FILES - - ${SOFASHELLS_SRC_DIR}/config.h.in - ${SOFASHELLS_SRC_DIR}/controller/MeshChangedEvent.h - ${SOFASHELLS_SRC_DIR}/controller/MeshInterpolator.h - ${SOFASHELLS_SRC_DIR}/controller/MeshInterpolator.inl - ${SOFASHELLS_SRC_DIR}/controller/TriangleSwitchExample.h - ${SOFASHELLS_SRC_DIR}/controller/TriangleSwitchExample.inl - ${SOFASHELLS_SRC_DIR}/engine/JoinMeshPoints.h - ${SOFASHELLS_SRC_DIR}/engine/JoinMeshPoints.inl - ${SOFASHELLS_SRC_DIR}/engine/FindClosePoints.h - ${SOFASHELLS_SRC_DIR}/engine/FindClosePoints.inl - ${SOFASHELLS_SRC_DIR}/forcefield/BezierTriangularBendingFEMForceField.h - ${SOFASHELLS_SRC_DIR}/forcefield/BezierTriangularBendingFEMForceField.inl - ${SOFASHELLS_SRC_DIR}/forcefield/CstFEMForceField.h - ${SOFASHELLS_SRC_DIR}/forcefield/CstFEMForceField.inl - ${SOFASHELLS_SRC_DIR}/forcefield/TriangularBendingFEMForceField.h - ${SOFASHELLS_SRC_DIR}/forcefield/TriangularBendingFEMForceField.inl - ${SOFASHELLS_SRC_DIR}/forcefield/TriangularShellForceField.h - ${SOFASHELLS_SRC_DIR}/forcefield/TriangularShellForceField.inl - ${SOFASHELLS_SRC_DIR}/mapping/BendingPlateMechanicalMapping.h - ${SOFASHELLS_SRC_DIR}/mapping/BendingPlateMechanicalMapping.inl - ${SOFASHELLS_SRC_DIR}/mapping/BezierTriangleMechanicalMapping.h - ${SOFASHELLS_SRC_DIR}/mapping/BezierTriangleMechanicalMapping.inl - ${SOFASHELLS_SRC_DIR}/misc/PointProjection.h - ${SOFASHELLS_SRC_DIR}/misc/PointProjection.inl - ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolation.h - ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolation.inl - ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolationM.h - ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolationM.inl - ${SOFASHELLS_SRC_DIR}/shells2/forcefield/BezierShellForceField.h - ${SOFASHELLS_SRC_DIR}/shells2/forcefield/BezierShellForceField.inl - ${SOFASHELLS_SRC_DIR}/shells2/mapping/BezierShellMechanicalMapping.h - ${SOFASHELLS_SRC_DIR}/shells2/mapping/BezierShellMechanicalMapping.inl - ) + ${SOFASHELLS_SRC_DIR}/config.h.in + ${SOFASHELLS_SRC_DIR}/controller/MeshChangedEvent.h + ${SOFASHELLS_SRC_DIR}/controller/MeshInterpolator.h + ${SOFASHELLS_SRC_DIR}/controller/MeshInterpolator.inl + ${SOFASHELLS_SRC_DIR}/controller/TriangleSwitchExample.h + ${SOFASHELLS_SRC_DIR}/controller/TriangleSwitchExample.inl + ${SOFASHELLS_SRC_DIR}/engine/JoinMeshPoints.h + ${SOFASHELLS_SRC_DIR}/engine/JoinMeshPoints.inl + ${SOFASHELLS_SRC_DIR}/engine/FindClosePoints.h + ${SOFASHELLS_SRC_DIR}/engine/FindClosePoints.inl + ${SOFASHELLS_SRC_DIR}/forcefield/BezierTriangularBendingFEMForceField.h + ${SOFASHELLS_SRC_DIR}/forcefield/BezierTriangularBendingFEMForceField.inl + ${SOFASHELLS_SRC_DIR}/forcefield/CstFEMForceField.h + ${SOFASHELLS_SRC_DIR}/forcefield/CstFEMForceField.inl + ${SOFASHELLS_SRC_DIR}/forcefield/TriangularBendingFEMForceField.h + ${SOFASHELLS_SRC_DIR}/forcefield/TriangularBendingFEMForceField.inl + ${SOFASHELLS_SRC_DIR}/forcefield/TriangularShellForceField.h + ${SOFASHELLS_SRC_DIR}/forcefield/TriangularShellForceField.inl + ${SOFASHELLS_SRC_DIR}/mapping/BendingPlateMechanicalMapping.h + ${SOFASHELLS_SRC_DIR}/mapping/BendingPlateMechanicalMapping.inl + ${SOFASHELLS_SRC_DIR}/mapping/BezierTriangleMechanicalMapping.h + ${SOFASHELLS_SRC_DIR}/mapping/BezierTriangleMechanicalMapping.inl + ${SOFASHELLS_SRC_DIR}/misc/PointProjection.h + ${SOFASHELLS_SRC_DIR}/misc/PointProjection.inl + ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolation.h + ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolation.inl + ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolationM.h + ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolationM.inl + ${SOFASHELLS_SRC_DIR}/shells2/forcefield/BezierShellForceField.h + ${SOFASHELLS_SRC_DIR}/shells2/forcefield/BezierShellForceField.inl + ${SOFASHELLS_SRC_DIR}/shells2/mapping/BezierShellMechanicalMapping.h + ${SOFASHELLS_SRC_DIR}/shells2/mapping/BezierShellMechanicalMapping.inl +) set(SOURCE_FILES - - ${SOFASHELLS_SRC_DIR}/initPluginShells.cpp - ${SOFASHELLS_SRC_DIR}/controller/MeshChangedEvent.cpp - ${SOFASHELLS_SRC_DIR}/controller/MeshInterpolator.cpp - ${SOFASHELLS_SRC_DIR}/controller/TriangleSwitchExample.cpp - ${SOFASHELLS_SRC_DIR}/engine/JoinMeshPoints.cpp - ${SOFASHELLS_SRC_DIR}/engine/FindClosePoints.cpp - ${SOFASHELLS_SRC_DIR}/forcefield/BezierTriangularBendingFEMForceField.cpp - ${SOFASHELLS_SRC_DIR}/forcefield/CstFEMForceField.cpp - ${SOFASHELLS_SRC_DIR}/forcefield/TriangularBendingFEMForceField.cpp - ${SOFASHELLS_SRC_DIR}/forcefield/TriangularShellForceField.cpp - ${SOFASHELLS_SRC_DIR}/mapping/BendingPlateMechanicalMapping.cpp - ${SOFASHELLS_SRC_DIR}/mapping/BezierTriangleMechanicalMapping.cpp - ${SOFASHELLS_SRC_DIR}/misc/PointProjection.cpp - ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolation.cpp - ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolationM.cpp - ${SOFASHELLS_SRC_DIR}/shells2/forcefield/BezierShellForceField.cpp - ${SOFASHELLS_SRC_DIR}/shells2/mapping/BezierShellMechanicalMapping.cpp - ) - - - - - + ${SOFASHELLS_SRC_DIR}/initPluginShells.cpp + ${SOFASHELLS_SRC_DIR}/controller/MeshChangedEvent.cpp + ${SOFASHELLS_SRC_DIR}/controller/MeshInterpolator.cpp + ${SOFASHELLS_SRC_DIR}/controller/TriangleSwitchExample.cpp + ${SOFASHELLS_SRC_DIR}/engine/JoinMeshPoints.cpp + ${SOFASHELLS_SRC_DIR}/engine/FindClosePoints.cpp + ${SOFASHELLS_SRC_DIR}/forcefield/BezierTriangularBendingFEMForceField.cpp + ${SOFASHELLS_SRC_DIR}/forcefield/CstFEMForceField.cpp + ${SOFASHELLS_SRC_DIR}/forcefield/TriangularBendingFEMForceField.cpp + ${SOFASHELLS_SRC_DIR}/forcefield/TriangularShellForceField.cpp + ${SOFASHELLS_SRC_DIR}/mapping/BendingPlateMechanicalMapping.cpp + ${SOFASHELLS_SRC_DIR}/mapping/BezierTriangleMechanicalMapping.cpp + ${SOFASHELLS_SRC_DIR}/misc/PointProjection.cpp + ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolation.cpp + ${SOFASHELLS_SRC_DIR}/shells2/fem/BezierShellInterpolationM.cpp + ${SOFASHELLS_SRC_DIR}/shells2/forcefield/BezierShellForceField.cpp + ${SOFASHELLS_SRC_DIR}/shells2/mapping/BezierShellMechanicalMapping.cpp +) if(SOFA-PLUGIN_SHELLS_ADAPTIVITY) - set(COMPILER_DEFINE "SOFA_BUILD_SHELLS_ADAPTIVITY") - - list(APPEND HEADER_FILES - - ${SOFASHELLS_SRC_DIR}/controller/AdaptiveCuttingController.h - ${SOFASHELLS_SRC_DIR}/controller/AdaptiveCuttingController.inl - ${SOFASHELLS_SRC_DIR}/controller/Test2DAdapter.h - ${SOFASHELLS_SRC_DIR}/controller/Test2DAdapter.inl - ${SOFASHELLS_SRC_DIR}/misc/Optimize2DSurface.h - ${SOFASHELLS_SRC_DIR}/misc/Optimize2DSurface.inl - ${SOFASHELLS_SRC_DIR}/misc/SurfaceParametrization.h - ${SOFASHELLS_SRC_DIR}/misc/SurfaceParametrization.inl - ) - - - list(APPEND SOURCE_FILES - - ${SOFASHELLS_SRC_DIR}/controller/AdaptiveCuttingController.cpp - ${SOFASHELLS_SRC_DIR}/controller/Test2DAdapter.cpp - ${SOFASHELLS_SRC_DIR}/misc/Optimize2DSurface.cpp - ${SOFASHELLS_SRC_DIR}/misc/SurfaceParametrization.cpp - ) - - - if(SofaGui_FOUND AND SofaOpenglVisual_FOUND) + set(COMPILER_DEFINE "SOFA_BUILD_SHELLS_ADAPTIVITY") + + list(APPEND HEADER_FILES + ${SOFASHELLS_SRC_DIR}/controller/AdaptiveCuttingController.h + ${SOFASHELLS_SRC_DIR}/controller/AdaptiveCuttingController.inl + ${SOFASHELLS_SRC_DIR}/controller/Test2DAdapter.h + ${SOFASHELLS_SRC_DIR}/controller/Test2DAdapter.inl + ${SOFASHELLS_SRC_DIR}/misc/Optimize2DSurface.h + ${SOFASHELLS_SRC_DIR}/misc/Optimize2DSurface.inl + ${SOFASHELLS_SRC_DIR}/misc/SurfaceParametrization.h + ${SOFASHELLS_SRC_DIR}/misc/SurfaceParametrization.inl + ) - list(APPEND HEADER_FILES - ${SOFASHELLS_SRC_DIR}/cutting/AdaptiveCutting.h - ) + list(APPEND SOURCE_FILES + ${SOFASHELLS_SRC_DIR}/controller/AdaptiveCuttingController.cpp + ${SOFASHELLS_SRC_DIR}/controller/Test2DAdapter.cpp + ${SOFASHELLS_SRC_DIR}/misc/Optimize2DSurface.cpp + ${SOFASHELLS_SRC_DIR}/misc/SurfaceParametrization.cpp + ) - list(APPEND SOURCE_FILES - ${SOFASHELLS_SRC_DIR}/cutting/AdaptiveCutting.cpp - ) + if(SofaGui_FOUND AND SofaOpenglVisual_FOUND) + list(APPEND HEADER_FILES + ${SOFASHELLS_SRC_DIR}/cutting/AdaptiveCutting.h + ) - endif() + list(APPEND SOURCE_FILES + ${SOFASHELLS_SRC_DIR}/cutting/AdaptiveCutting.cpp + ) + endif() endif() @@ -131,21 +116,27 @@ endif() # Create the plugin library add_library(${PROJECT_NAME} SHARED ${HEADER_FILES} ${SOURCE_FILES} ${README_FILES}) -# Link the plugin library to its dependency -target_link_libraries(${PROJECT_NAME} Sofa.Component.Controller Sofa.Helper Sofa.Type Sofa.Component.Topology.Container.Dynamic Sofa.Component.StateContainer - Sofa.Core Sofa.Geometry Sofa.GL Sofa.Component.Collision.Detection.Intersection Sofa.Component.Mapping.NonLinear Sofa.Component.Constraint.Lagrangian.Model - Sofa.Component.Constraint.Lagrangian.Solver Sofa.Component.Constraint.Lagrangian Sofa.Component.Constraint - Sofa.Simulation.Core Sofa.Component.Constraint.Projective Sofa.Component.Mass Sofa.Component.SolidMechanics.Spring Sofa.Component.MechanicalLoad - Sofa.Component.LinearSolver.Iterative Sofa.Component.ODESolver.Backward Sofa.Component.Engine.Select Sofa.Component.Mapping.Linear +# Link the plugin library to its dependency(ies). +target_link_libraries(${PROJECT_NAME} + Sofa.Component.Controller + Sofa.Component.Topology.Container.Dynamic + Sofa.Component.StateContainer + Sofa.GL + Sofa.Component.Collision.Detection.Intersection + Sofa.Component.Mapping.NonLinear + Sofa.Component.Constraint.Lagrangian.Model + Sofa.Component.Constraint.Lagrangian.Solver + Sofa.Component.Constraint.Lagrangian + Sofa.Component.Constraint.Projective + Sofa.Component.Mass + Sofa.Component.SolidMechanics.Spring + Sofa.Component.MechanicalLoad + Sofa.Component.LinearSolver.Iterative + Sofa.Component.ODESolver.Backward + Sofa.Component.Engine.Select + Sofa.Component.Mapping.Linear ) - -if(SofaGui_FOUND AND SofaOpenglVisual_FOUND) - target_link_libraries(${PROJECT_NAME} SofaGui SofaOpenglVisual) -endif() - - -# Create package Config, Version & Target files sofa_create_package_with_targets( PACKAGE_NAME ${PROJECT_NAME} PACKAGE_VERSION ${PROJECT_VERSION} diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..f59945c --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,504 @@ + GNU LESSER GENERAL PUBLIC LICENSE + Version 2.1, February 1999 + + Copyright (C) 1991, 1999 Free Software Foundation, Inc. + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + +[This is the first released version of the Lesser GPL. It also counts + as the successor of the GNU Library Public License, version 2, hence + the version number 2.1.] + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +Licenses are intended to guarantee your freedom to share and change +free software--to make sure the software is free for all its users. + + This license, the Lesser General Public License, applies to some +specially designated software packages--typically libraries--of the +Free Software Foundation and other authors who decide to use it. You +can use it too, but we suggest you first think carefully about whether +this license or the ordinary General Public License is the better +strategy to use in any particular case, based on the explanations below. + + When we speak of free software, we are referring to freedom of use, +not price. Our General Public Licenses are designed to make sure that +you have the freedom to distribute copies of free software (and charge +for this service if you wish); that you receive source code or can get +it if you want it; that you can change the software and use pieces of +it in new free programs; and that you are informed that you can do +these things. + + To protect your rights, we need to make restrictions that forbid +distributors to deny you these rights or to ask you to surrender these +rights. These restrictions translate to certain responsibilities for +you if you distribute copies of the library or if you modify it. + + For example, if you distribute copies of the library, whether gratis +or for a fee, you must give the recipients all the rights that we gave +you. You must make sure that they, too, receive or can get the source +code. If you link other code with the library, you must provide +complete object files to the recipients, so that they can relink them +with the library after making changes to the library and recompiling +it. And you must show them these terms so they know their rights. + + We protect your rights with a two-step method: (1) we copyright the +library, and (2) we offer you this license, which gives you legal +permission to copy, distribute and/or modify the library. + + To protect each distributor, we want to make it very clear that +there is no warranty for the free library. Also, if the library is +modified by someone else and passed on, the recipients should know +that what they have is not the original version, so that the original +author's reputation will not be affected by problems that might be +introduced by others. + + Finally, software patents pose a constant threat to the existence of +any free program. We wish to make sure that a company cannot +effectively restrict the users of a free program by obtaining a +restrictive license from a patent holder. Therefore, we insist that +any patent license obtained for a version of the library must be +consistent with the full freedom of use specified in this license. + + Most GNU software, including some libraries, is covered by the +ordinary GNU General Public License. This license, the GNU Lesser +General Public License, applies to certain designated libraries, and +is quite different from the ordinary General Public License. We use +this license for certain libraries in order to permit linking those +libraries into non-free programs. + + When a program is linked with a library, whether statically or using +a shared library, the combination of the two is legally speaking a +combined work, a derivative of the original library. The ordinary +General Public License therefore permits such linking only if the +entire combination fits its criteria of freedom. The Lesser General +Public License permits more lax criteria for linking other code with +the library. + + We call this license the "Lesser" General Public License because it +does Less to protect the user's freedom than the ordinary General +Public License. It also provides other free software developers Less +of an advantage over competing non-free programs. These disadvantages +are the reason we use the ordinary General Public License for many +libraries. However, the Lesser license provides advantages in certain +special circumstances. + + For example, on rare occasions, there may be a special need to +encourage the widest possible use of a certain library, so that it becomes +a de-facto standard. To achieve this, non-free programs must be +allowed to use the library. A more frequent case is that a free +library does the same job as widely used non-free libraries. In this +case, there is little to gain by limiting the free library to free +software only, so we use the Lesser General Public License. + + In other cases, permission to use a particular library in non-free +programs enables a greater number of people to use a large body of +free software. For example, permission to use the GNU C Library in +non-free programs enables many more people to use the whole GNU +operating system, as well as its variant, the GNU/Linux operating +system. + + Although the Lesser General Public License is Less protective of the +users' freedom, it does ensure that the user of a program that is +linked with the Library has the freedom and the wherewithal to run +that program using a modified version of the Library. + + The precise terms and conditions for copying, distribution and +modification follow. Pay close attention to the difference between a +"work based on the library" and a "work that uses the library". The +former contains code derived from the library, whereas the latter must +be combined with the library in order to run. + + GNU LESSER GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License Agreement applies to any software library or other +program which contains a notice placed by the copyright holder or +other authorized party saying it may be distributed under the terms of +this Lesser General Public License (also called "this License"). +Each licensee is addressed as "you". + + A "library" means a collection of software functions and/or data +prepared so as to be conveniently linked with application programs +(which use some of those functions and data) to form executables. + + The "Library", below, refers to any such software library or work +which has been distributed under these terms. A "work based on the +Library" means either the Library or any derivative work under +copyright law: that is to say, a work containing the Library or a +portion of it, either verbatim or with modifications and/or translated +straightforwardly into another language. (Hereinafter, translation is +included without limitation in the term "modification".) + + "Source code" for a work means the preferred form of the work for +making modifications to it. For a library, complete source code means +all the source code for all modules it contains, plus any associated +interface definition files, plus the scripts used to control compilation +and installation of the library. + + Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running a program using the Library is not restricted, and output from +such a program is covered only if its contents constitute a work based +on the Library (independent of the use of the Library in a tool for +writing it). Whether that is true depends on what the Library does +and what the program that uses the Library does. + + 1. You may copy and distribute verbatim copies of the Library's +complete source code as you receive it, in any medium, provided that +you conspicuously and appropriately publish on each copy an +appropriate copyright notice and disclaimer of warranty; keep intact +all the notices that refer to this License and to the absence of any +warranty; and distribute a copy of this License along with the +Library. + + You may charge a fee for the physical act of transferring a copy, +and you may at your option offer warranty protection in exchange for a +fee. + + 2. You may modify your copy or copies of the Library or any portion +of it, thus forming a work based on the Library, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) The modified work must itself be a software library. + + b) You must cause the files modified to carry prominent notices + stating that you changed the files and the date of any change. + + c) You must cause the whole of the work to be licensed at no + charge to all third parties under the terms of this License. + + d) If a facility in the modified Library refers to a function or a + table of data to be supplied by an application program that uses + the facility, other than as an argument passed when the facility + is invoked, then you must make a good faith effort to ensure that, + in the event an application does not supply such function or + table, the facility still operates, and performs whatever part of + its purpose remains meaningful. + + (For example, a function in a library to compute square roots has + a purpose that is entirely well-defined independent of the + application. Therefore, Subsection 2d requires that any + application-supplied function or table used by this function must + be optional: if the application does not supply it, the square + root function must still compute square roots.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Library, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Library, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote +it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Library. + +In addition, mere aggregation of another work not based on the Library +with the Library (or with a work based on the Library) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may opt to apply the terms of the ordinary GNU General Public +License instead of this License to a given copy of the Library. To do +this, you must alter all the notices that refer to this License, so +that they refer to the ordinary GNU General Public License, version 2, +instead of to this License. (If a newer version than version 2 of the +ordinary GNU General Public License has appeared, then you can specify +that version instead if you wish.) Do not make any other change in +these notices. + + Once this change is made in a given copy, it is irreversible for +that copy, so the ordinary GNU General Public License applies to all +subsequent copies and derivative works made from that copy. + + This option is useful when you wish to copy part of the code of +the Library into a program that is not a library. + + 4. You may copy and distribute the Library (or a portion or +derivative of it, under Section 2) in object code or executable form +under the terms of Sections 1 and 2 above provided that you accompany +it with the complete corresponding machine-readable source code, which +must be distributed under the terms of Sections 1 and 2 above on a +medium customarily used for software interchange. + + If distribution of object code is made by offering access to copy +from a designated place, then offering equivalent access to copy the +source code from the same place satisfies the requirement to +distribute the source code, even though third parties are not +compelled to copy the source along with the object code. + + 5. A program that contains no derivative of any portion of the +Library, but is designed to work with the Library by being compiled or +linked with it, is called a "work that uses the Library". Such a +work, in isolation, is not a derivative work of the Library, and +therefore falls outside the scope of this License. + + However, linking a "work that uses the Library" with the Library +creates an executable that is a derivative of the Library (because it +contains portions of the Library), rather than a "work that uses the +library". The executable is therefore covered by this License. +Section 6 states terms for distribution of such executables. + + When a "work that uses the Library" uses material from a header file +that is part of the Library, the object code for the work may be a +derivative work of the Library even though the source code is not. +Whether this is true is especially significant if the work can be +linked without the Library, or if the work is itself a library. The +threshold for this to be true is not precisely defined by law. + + If such an object file uses only numerical parameters, data +structure layouts and accessors, and small macros and small inline +functions (ten lines or less in length), then the use of the object +file is unrestricted, regardless of whether it is legally a derivative +work. (Executables containing this object code plus portions of the +Library will still fall under Section 6.) + + Otherwise, if the work is a derivative of the Library, you may +distribute the object code for the work under the terms of Section 6. +Any executables containing that work also fall under Section 6, +whether or not they are linked directly with the Library itself. + + 6. As an exception to the Sections above, you may also combine or +link a "work that uses the Library" with the Library to produce a +work containing portions of the Library, and distribute that work +under terms of your choice, provided that the terms permit +modification of the work for the customer's own use and reverse +engineering for debugging such modifications. + + You must give prominent notice with each copy of the work that the +Library is used in it and that the Library and its use are covered by +this License. You must supply a copy of this License. If the work +during execution displays copyright notices, you must include the +copyright notice for the Library among them, as well as a reference +directing the user to the copy of this License. Also, you must do one +of these things: + + a) Accompany the work with the complete corresponding + machine-readable source code for the Library including whatever + changes were used in the work (which must be distributed under + Sections 1 and 2 above); and, if the work is an executable linked + with the Library, with the complete machine-readable "work that + uses the Library", as object code and/or source code, so that the + user can modify the Library and then relink to produce a modified + executable containing the modified Library. (It is understood + that the user who changes the contents of definitions files in the + Library will not necessarily be able to recompile the application + to use the modified definitions.) + + b) Use a suitable shared library mechanism for linking with the + Library. A suitable mechanism is one that (1) uses at run time a + copy of the library already present on the user's computer system, + rather than copying library functions into the executable, and (2) + will operate properly with a modified version of the library, if + the user installs one, as long as the modified version is + interface-compatible with the version that the work was made with. + + c) Accompany the work with a written offer, valid for at + least three years, to give the same user the materials + specified in Subsection 6a, above, for a charge no more + than the cost of performing this distribution. + + d) If distribution of the work is made by offering access to copy + from a designated place, offer equivalent access to copy the above + specified materials from the same place. + + e) Verify that the user has already received a copy of these + materials or that you have already sent this user a copy. + + For an executable, the required form of the "work that uses the +Library" must include any data and utility programs needed for +reproducing the executable from it. However, as a special exception, +the materials to be distributed need not include anything that is +normally distributed (in either source or binary form) with the major +components (compiler, kernel, and so on) of the operating system on +which the executable runs, unless that component itself accompanies +the executable. + + It may happen that this requirement contradicts the license +restrictions of other proprietary libraries that do not normally +accompany the operating system. Such a contradiction means you cannot +use both them and the Library together in an executable that you +distribute. + + 7. You may place library facilities that are a work based on the +Library side-by-side in a single library together with other library +facilities not covered by this License, and distribute such a combined +library, provided that the separate distribution of the work based on +the Library and of the other library facilities is otherwise +permitted, and provided that you do these two things: + + a) Accompany the combined library with a copy of the same work + based on the Library, uncombined with any other library + facilities. This must be distributed under the terms of the + Sections above. + + b) Give prominent notice with the combined library of the fact + that part of it is a work based on the Library, and explaining + where to find the accompanying uncombined form of the same work. + + 8. You may not copy, modify, sublicense, link with, or distribute +the Library except as expressly provided under this License. Any +attempt otherwise to copy, modify, sublicense, link with, or +distribute the Library is void, and will automatically terminate your +rights under this License. However, parties who have received copies, +or rights, from you under this License will not have their licenses +terminated so long as such parties remain in full compliance. + + 9. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Library or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Library (or any work based on the +Library), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Library or works based on it. + + 10. Each time you redistribute the Library (or any work based on the +Library), the recipient automatically receives a license from the +original licensor to copy, distribute, link with or modify the Library +subject to these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties with +this License. + + 11. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Library at all. For example, if a patent +license would not permit royalty-free redistribution of the Library by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Library. + +If any portion of this section is held invalid or unenforceable under any +particular circumstance, the balance of the section is intended to apply, +and the section as a whole is intended to apply in other circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 12. If the distribution and/or use of the Library is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Library under this License may add +an explicit geographical distribution limitation excluding those countries, +so that distribution is permitted only in or among countries not thus +excluded. In such case, this License incorporates the limitation as if +written in the body of this License. + + 13. The Free Software Foundation may publish revised and/or new +versions of the Lesser General Public License from time to time. +Such new versions will be similar in spirit to the present version, +but may differ in detail to address new problems or concerns. + +Each version is given a distinguishing version number. If the Library +specifies a version number of this License which applies to it and +"any later version", you have the option of following the terms and +conditions either of that version or of any later version published by +the Free Software Foundation. If the Library does not specify a +license version number, you may choose any version ever published by +the Free Software Foundation. + + 14. If you wish to incorporate parts of the Library into other free +programs whose distribution conditions are incompatible with these, +write to the author to ask for permission. For software which is +copyrighted by the Free Software Foundation, write to the Free +Software Foundation; we sometimes make exceptions for this. Our +decision will be guided by the two goals of preserving the free status +of all derivatives of our free software and of promoting the sharing +and reuse of software generally. + + NO WARRANTY + + 15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO +WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW. +EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR +OTHER PARTIES PROVIDE THE LIBRARY "AS IS" WITHOUT WARRANTY OF ANY +KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE +LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME +THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN +WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY +AND/OR REDISTRIBUTE THE LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU +FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR +CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE +LIBRARY (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING +RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A +FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF +SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH +DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Libraries + + If you develop a new library, and you want it to be of the greatest +possible use to the public, we recommend making it free software that +everyone can redistribute and change. You can do so by permitting +redistribution under these terms (or, alternatively, under the terms of the +ordinary General Public License). + + To apply these terms, attach the following notices to the library. It is +safest to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least the +"copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 + USA + +Also add information on how to contact you by electronic and paper mail. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the library, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the + library `Frob' (a library for tweaking knobs) written by James Random + Hacker. + + , 1 April 1990 + Ty Coon, President of Vice + +That's all there is to it! diff --git a/README.md b/README.md new file mode 100644 index 0000000..8b5f84a --- /dev/null +++ b/README.md @@ -0,0 +1,16 @@ +# Shell + +This plugin for [SOFA](https://github.com/sofa-framework/sofa) implements mechanical models following the Shell method. + +![Shell example](doc/shell.png) + +## AUTHORS + - Thomas Golembiovsky + - Christian Duriez + - Olivier Comas + - Igor Peterlik + - Stéphane Cotin + + +## LICENCE + - LGPL diff --git a/Shell.dox b/Shell.dox new file mode 100644 index 0000000..7055814 --- /dev/null +++ b/Shell.dox @@ -0,0 +1,5 @@ +/** @mainpage + +This plugin for SOFA implements mechanical models following the Shell method. + +*/ diff --git a/ShellConfig.cmake.in b/ShellConfig.cmake.in new file mode 100644 index 0000000..63e71a6 --- /dev/null +++ b/ShellConfig.cmake.in @@ -0,0 +1,13 @@ +# CMake package configuration file for the plugin @PROJECT_NAME@ + +@PACKAGE_GUARD@ +@PACKAGE_INIT@ + +find_package(Sofa.Component.Controller QUIET REQUIRED) +find_package(Sofa.Component.StateContainer QUIET REQUIRED) + +if(NOT TARGET @PROJECT_NAME@) + include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@Targets.cmake") +endif() + +check_required_components(@PROJECT_NAME@) diff --git a/doc/index.html b/doc/index.html new file mode 100644 index 0000000..e69de29 diff --git a/doc/shell.png b/doc/shell.png new file mode 100644 index 0000000..f3a0004 Binary files /dev/null and b/doc/shell.png differ diff --git a/examples/sofapython3/TriangularBendingFEMForceField.py b/examples/sofapython3/TriangularBendingFEMForceField.py index ca79297..975aaf1 100644 --- a/examples/sofapython3/TriangularBendingFEMForceField.py +++ b/examples/sofapython3/TriangularBendingFEMForceField.py @@ -1,14 +1,13 @@ def createScene(rootNode): + rootNode.addObject('RequiredPlugin', name="Shell") rootNode.gravity = [0, -98.1, 0] rootNode.addObject('AttachBodyButtonSetting', stiffness=0.1) - rootNode.addObject('FreeMotionAnimationLoop') - rootNode.addObject('GenericConstraintSolver', maxIterations=1e3, tolerance=1e-3) + rootNode.addObject('DefaultAnimationLoop') square = rootNode.addChild('Square') square.addObject('EulerImplicitSolver') square.addObject('SparseLDLSolver') - square.addObject('GenericConstraintCorrection') square.addObject('MeshOBJLoader', filename='mesh/square1.obj') square.addObject('MeshTopology', src=square.MeshOBJLoader.getLinkPath()) square.addObject('MechanicalObject', template='Rigid3') diff --git a/examples/xml/ShellTest.html b/examples/xml/ShellTest.html new file mode 100644 index 0000000..f7df77d --- /dev/null +++ b/examples/xml/ShellTest.html @@ -0,0 +1,13 @@ + + + + + +Shell +

Shell

+
This example shows a simple simulation, as shown in the following image.
+More documentation can be found in doc/index.html
+
+
+
+ diff --git a/examples/xml/ShellTest.scn b/examples/xml/ShellTest.scn new file mode 100644 index 0000000..95a6ac0 --- /dev/null +++ b/examples/xml/ShellTest.scn @@ -0,0 +1,45 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/Shell/config.h.in b/src/Shell/config.h.in new file mode 100644 index 0000000..b0e62d0 --- /dev/null +++ b/src/Shell/config.h.in @@ -0,0 +1,34 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +#define SHELL_VERSION @PROJECT_VERSION@ + +#ifdef SOFA_BUILD_SHELL +# define SOFA_TARGET @PROJECT_NAME@ +# define SOFA_SHELL_API SOFA_EXPORT_DYNAMIC_LIBRARY +#else +# define SOFA_SHELL_API SOFA_IMPORT_DYNAMIC_LIBRARY +#endif + diff --git a/src/Shell/controller/MeshChangedEvent.cpp b/src/Shell/controller/MeshChangedEvent.cpp new file mode 100644 index 0000000..9138366 --- /dev/null +++ b/src/Shell/controller/MeshChangedEvent.cpp @@ -0,0 +1,30 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ + +#include + +namespace shell::objectmodel +{ + +SOFA_EVENT_CPP( MeshChangedEvent ) + +} // namespace diff --git a/src/Shell/controller/MeshChangedEvent.h b/src/Shell/controller/MeshChangedEvent.h new file mode 100644 index 0000000..ffdcf38 --- /dev/null +++ b/src/Shell/controller/MeshChangedEvent.h @@ -0,0 +1,53 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include + +namespace shell::objectmodel +{ + +using namespace sofa::type; + +// Sent by MeshInterpolator when the mesh changes +class MeshChangedEvent : public sofa::core::objectmodel::Event +{ +public: + + SOFA_EVENT_H( MeshChangedEvent ) + + MeshChangedEvent(double _alpha) : alpha(_alpha) + {} + + ~MeshChangedEvent() {} + + double getAlpha() const { return alpha; } + +private: + + double alpha; // interpolation parameter +}; + +} // namespace diff --git a/src/Shell/controller/MeshInterpolator.cpp b/src/Shell/controller/MeshInterpolator.cpp new file mode 100644 index 0000000..154eac4 --- /dev/null +++ b/src/Shell/controller/MeshInterpolator.cpp @@ -0,0 +1,39 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SHELL_CONTROLLER_MESHINTERPOLATION_CPP + +#include +#include + +namespace shell::controller +{ + +using namespace sofa::defaulttype; + +// Register in the Factory +int MeshInterpolatorClass = sofa::core::RegisterObject("Performs linear interpolation between two meshes") +.add< MeshInterpolator >(true) // default template +; + +template class SOFA_SHELL_API MeshInterpolator; + +} // namespace diff --git a/src/Shell/controller/MeshInterpolator.h b/src/Shell/controller/MeshInterpolator.h new file mode 100644 index 0000000..14748fb --- /dev/null +++ b/src/Shell/controller/MeshInterpolator.h @@ -0,0 +1,82 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include + +#include + +namespace shell::controller +{ + +template +class MeshInterpolator : public sofa::component::controller::Controller +{ +public: + SOFA_CLASS(SOFA_TEMPLATE(MeshInterpolator,DataTypes), sofa::component::controller::Controller); + + typedef typename DataTypes::VecCoord VecCoord; + typedef typename DataTypes::VecDeriv VecDeriv; + typedef typename DataTypes::Coord Coord; + typedef typename DataTypes::Deriv Deriv; + typedef typename Coord::value_type Real; + typedef sofa::type::Vec<3,Real> Vec3; + +protected: + + MeshInterpolator(); + + virtual ~MeshInterpolator(); + +public: + + void init() override; + void reinit() override; + + sofa::Data f_startTime; + sofa::Data f_nbSteps; + sofa::Data f_increment; + + sofa::Data f_startPosition; + sofa::Data > f_startNormals; + sofa::Data f_endPosition; + sofa::Data > f_endNormals; + + sofa::Data f_position; + sofa::Data > f_normals; + + void onEndAnimationStep(const double dt) override; + + Real getInterpolationVar() { return alpha; } + +private: + + unsigned int stepCounter; + Real alpha; + + void interpolate(); +}; + +} // namespace diff --git a/src/Shell/controller/MeshInterpolator.inl b/src/Shell/controller/MeshInterpolator.inl new file mode 100644 index 0000000..ea02844 --- /dev/null +++ b/src/Shell/controller/MeshInterpolator.inl @@ -0,0 +1,191 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include + + +namespace shell::controller +{ + +template +MeshInterpolator::MeshInterpolator() +: f_startTime(initData(&f_startTime, (Real)0, "startTime", "Time when the interpolation starts")) +, f_nbSteps(initData(&f_nbSteps, (unsigned int)10, "nbSteps", "Perform a transition every nbStep steps")) +, f_increment(initData(&f_increment, (Real)0.05, "increment", "How quickly converge to the final positions")) +, f_startPosition(initData(&f_startPosition, "startPosition","Starting positions of the nodes")) +, f_startNormals(initData(&f_startNormals, "startNormals","Starting normals of the nodes")) +, f_endPosition(initData(&f_endPosition, "endPosition","Final positions of the nodes")) +, f_endNormals(initData(&f_endNormals, "endNormals","Final normals of the nodes")) +, f_position(initData(&f_position, "position","Interpolated positions of the nodes")) +, f_normals(initData(&f_normals, "normals","Interpolated normals of the nodes")) +, stepCounter(0) +{ +} + +template +MeshInterpolator::~MeshInterpolator() +{ +} + +template +void MeshInterpolator::init() +{ + reinit(); +} + +template +void MeshInterpolator::reinit() +{ + *this->f_listening.beginEdit() = true; + this->f_listening.endEdit(); + + unsigned int lenStart = f_startPosition.getValue().size(); + unsigned int lenEnd = f_endPosition.getValue().size(); + + // Check that the number of nodes is the same + if (lenStart != lenEnd) { + msg_warning() << "Number of start and end nodes has to be the same."; + if (lenStart > lenEnd) { + VecCoord &pts = *f_startPosition.beginEdit(); + pts.resize(lenEnd); + f_startPosition.endEdit(); + lenStart = lenEnd; + } else { + VecCoord &pts = *f_endPosition.beginEdit(); + pts.resize(lenStart); + f_endPosition.endEdit(); + } + } + + unsigned int lenStartN = f_startNormals.getValue().size(); + unsigned int lenEndN = f_endNormals.getValue().size(); + + // Check that the number of nodes is the same + if (((lenStartN != 0) || (lenEndN != 0)) && + ((lenStartN != lenEndN) || (lenStart != lenStartN))) { + msg_warning() << "Number of start normals, end normals and positions has to be the same."; + f_startNormals.beginEdit()->clear(); + f_startNormals.endEdit(); + f_endNormals.beginEdit()->clear(); + f_endNormals.endEdit(); + } + + // Check startTime + if (f_startTime.getValue() < 0.0) { + msg_warning() << "startTime has to be greater then or equal to 0."; + *f_startTime.beginEdit() = 0.0; + f_startTime.endEdit(); + } + + // Check nbSteps + if (f_nbSteps.getValue() == 0) { + msg_warning() << "nbSteps has to be nonzero."; + *f_nbSteps.beginEdit() = 1; + f_nbSteps.endEdit(); + } + + if ((f_increment.getValue() <= 0.0) || (f_increment.getValue() > 1.0)) { + msg_warning() << "Increment has to be geater than 0 and " + "smaller than or equall to 1."; + *f_increment.beginEdit() = 0.05; + f_increment.endEdit(); + } + + // XXX: does it make sense to reinit also internal state? + stepCounter = 0; + alpha = 0; + + // Start with starting point + *f_position.beginEdit() = f_startPosition.getValue(); + f_position.endEdit(); + + // Start with starting normals + *f_normals.beginEdit() = f_startNormals.getValue(); + f_normals.endEdit(); + +} + + +template +void MeshInterpolator::onEndAnimationStep(const double /*dt*/) +{ + if (alpha >= 1.0) + return; // Nothing more to do + + if (getContext()->getTime() < f_startTime.getValue()) + return; // Not yet ... + + stepCounter++; + + if (stepCounter < f_nbSteps.getValue()) + return; // Stay still for a few steps + + stepCounter = 0; + + // Increase the linear factor + alpha += f_increment.getValue(); + + // Update positions + interpolate(); + + // Send signal about change + shell::objectmodel::MeshChangedEvent mcEvent(alpha); + this->getContext()->propagateEvent(sofa::core::execparams::defaultInstance(), &mcEvent); +} + +template +void MeshInterpolator::interpolate() +{ + const VecCoord &startPt = f_startPosition.getValue(); + const VecCoord &endPt = f_endPosition.getValue(); + + const sofa::type::vector &startNorm = f_startNormals.getValue(); + const sofa::type::vector &endNorm = f_endNormals.getValue(); + + VecCoord &pt = *f_position.beginEdit(); + sofa::type::vector &norm = *f_normals.beginEdit(); + + for (unsigned int i=0; i 0) { + for (unsigned int i=0; i. * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SHELL_ENGINE_JOINMESHPOINTS_CPP + +#include +#include + +namespace shell::engine +{ + +int JoinMeshPointsClass = sofa::core::RegisterObject( + "Join two or more points of the mesh together while updating the topology" + " accordingly") + +.add< JoinMeshPoints >(true) // default template + +.add< JoinMeshPoints >() +.add< JoinMeshPoints >() +.add< JoinMeshPoints >() +.add< JoinMeshPoints >() +; + +template class SOFA_SHELL_API JoinMeshPoints; +template class SOFA_SHELL_API JoinMeshPoints; +template class SOFA_SHELL_API JoinMeshPoints; +template class SOFA_SHELL_API JoinMeshPoints; +template class SOFA_SHELL_API JoinMeshPoints; + +} // namespace diff --git a/src/Shell/engine/JoinMeshPoints.h b/src/Shell/engine/JoinMeshPoints.h new file mode 100644 index 0000000..114bd4f --- /dev/null +++ b/src/Shell/engine/JoinMeshPoints.h @@ -0,0 +1,132 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace shell::engine +{ + +/** + * @brief Joins two or more points of the mesh together and updates the + * topology accordingly. + * + * @tparam DataTypes Associated data type. + */ +template +class JoinMeshPoints : public sofa::core::DataEngine +{ +public: + SOFA_CLASS(SOFA_TEMPLATE(JoinMeshPoints,DataTypes), sofa::core::DataEngine); + + typedef typename DataTypes::VecCoord VecCoord; + typedef typename DataTypes::Coord Coord; + typedef typename Coord::value_type Real; + + typedef sofa::type::Vec<3,Real> Vec3; + typedef unsigned int Index; + +protected: + JoinMeshPoints(); + + virtual ~JoinMeshPoints(); + +public: + + void init() override; + void reinit() override; + void doUpdate() override; + + // TODO: add methods to find out the inverse mappings + + /** + * @brief Get index of a node in input topology based on the triangle and + * node ID in output topology. + * + * @param triangleId Triangle from output topology. + * @param nodeId Node in output topology. + * + * @return Node from input topology. + */ + Index getSrcNodeFromTri(Index triangleId, Index nodeId) { + if (triangleId < f_output_triangles.getValue().size()) { + if (nodeId < f_output_position.getValue().size()) { + const std::map &n2n = imapTriNode2Node[triangleId]; + std::map::const_iterator iter = n2n.find(nodeId); + if (iter != n2n.end()) { + return iter->second; + } else { + // Not a joined point + return imapNode2Node[nodeId]; + } + } + } + // Invalid indices + msg_error() << "Requested invalid triangle " << triangleId << " and node " << nodeId; + return 0; + } + + + // Data + + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_input_joinPoints; + + sofa::Data f_input_position; + sofa::Data< sofa::type::vector > f_input_normals; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_input_edges; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_input_triangles; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_input_quads; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_input_tetrahedra; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_input_hexahedra; + + sofa::Data f_output_position; + sofa::Data< sofa::type::vector > f_output_normals; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_output_edges; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_output_triangles; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_output_quads; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_output_tetrahedra; + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > f_output_hexahedra; + sofa::Data f_output_mergedPosition; + sofa::Data< sofa::type::vector > f_output_mergedNormals; + +private: + + template void createElements( + std::map mapInIn, + sofa::type::vector mapInOut, + const sofa::Data< sofa::type::vector< sofa::type::fixed_array > > &inElements, + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > &outElements); + + // Inverse mappings: from output to input topology + // ... node -> node (only for nodes that are not joined and thus have one to one mapping) + std::map imapNode2Node; + // ... triangle + node -> node + std::map > imapTriNode2Node; +}; + +} // namespace diff --git a/src/Shell/engine/JoinMeshPoints.inl b/src/Shell/engine/JoinMeshPoints.inl new file mode 100644 index 0000000..e0611df --- /dev/null +++ b/src/Shell/engine/JoinMeshPoints.inl @@ -0,0 +1,272 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include + +namespace shell::engine +{ + +using namespace sofa::type; +using namespace sofa::core::objectmodel; + +template +JoinMeshPoints::JoinMeshPoints() +: f_input_joinPoints(initData(&f_input_joinPoints,"joinPoints","Which points to join")) + , f_input_position(initData(&f_input_position,"position","Input Vertices")) + , f_input_normals(initData(&f_input_normals,"normals","Input Normals")) + , f_input_edges(initData(&f_input_edges,"edges","Input Edges")) + , f_input_triangles(initData(&f_input_triangles,"triangles","Input Triangles")) + , f_input_quads(initData(&f_input_quads,"quads","Input Quads")) + , f_input_tetrahedra(initData(&f_input_tetrahedra,"tetrahedra","Input Tetrahedra")) + , f_input_hexahedra(initData(&f_input_hexahedra,"hexahedra","Input Hexahedra")) + + , f_output_position(initData(&f_output_position,"joinedPosition","Output Vertices of the joined mesh")) + , f_output_normals(initData(&f_output_normals,"joinedNormals","Output Normals of the joined mesh")) + , f_output_edges(initData(&f_output_edges,"joinedEdges","Output Edges of the joined mesh")) + , f_output_triangles(initData(&f_output_triangles,"joinedTriangles","Output Triangles of the joined mesh")) + , f_output_quads(initData(&f_output_quads,"joinedQuads","Output Quads of the joined mesh")) + , f_output_tetrahedra(initData(&f_output_tetrahedra,"joinTetrahedra","Output Tetrahedra of the joined mesh")) + , f_output_hexahedra(initData(&f_output_hexahedra,"joinHexahedra","Output Hexahedra of the joined mesh")) + , f_output_mergedPosition(initData(&f_output_mergedPosition,"mergedPosition","Positions of the merged Vertices of the input topology")) + , f_output_mergedNormals(initData(&f_output_mergedNormals,"mergedNormals","Normals of the merged Vertices of the input topology")) +{ +} + +template +JoinMeshPoints::~JoinMeshPoints() +{ +} + +template +void JoinMeshPoints::init() +{ + addInput(&f_input_joinPoints); + + addInput(&f_input_position); + addInput(&f_input_normals); + addInput(&f_input_edges); + addInput(&f_input_triangles); + addInput(&f_input_quads); + addInput(&f_input_tetrahedra); + addInput(&f_input_hexahedra); + + addOutput(&f_output_position); + addOutput(&f_output_normals); + addOutput(&f_output_edges); + addOutput(&f_output_triangles); + addOutput(&f_output_quads); + addOutput(&f_output_tetrahedra); + addOutput(&f_output_hexahedra); + + addOutput(&f_output_mergedPosition); + addOutput(&f_output_mergedNormals); + + setDirtyValue(); +} + +template +void JoinMeshPoints::reinit() +{ + update(); +} + +template +void JoinMeshPoints::doUpdate() +{ + + const sofa::type::vector< sofa::type::fixed_array >& inJP = f_input_joinPoints.getValue(); + const VecCoord& inPt = f_input_position.getValue(); + const sofa::type::vector& inNorm = f_input_normals.getValue(); + + VecCoord& outPt = *f_output_position.beginEdit(); + sofa::type::vector& outNorm = *f_output_normals.beginEdit(); + VecCoord& outMPt = *f_output_mergedPosition.beginEdit(); + sofa::type::vector& outMNorm = *f_output_mergedNormals.beginEdit(); + outPt.clear(); + outNorm.clear(); + outMPt.resize(inPt.size()); + outMNorm.resize(inPt.size()); + + bool bHasNormals = false; + if (inPt.size() == inNorm.size()) { + bHasNormals = true; + } else if (inNorm.size() != 0) { + msg_warning() << "Normal count does not match node count! Ignoring normals."; + } + + // Map analogy to the value of f_input_joinPoints. It maps index of a node + // from input array to index into output array. + std::map mapInIn; + + // Add items so that if i maps to j then i > j + for (Index i=0; i::iterator iter = + mapInIn.find(from); + if (iter != mapInIn.end()) { + if (iter->second == to) { + continue; + } else if (iter->second < to) { + // Reassign the current item and keep the old one untact + from = to; + to = iter->second; + } else { + // Reassign previous item + Index newFrom = iter->second; + Index newTo = to; + mapInIn.erase(iter); + mapInIn[newFrom] = newTo; + // ... and add the new item + } + } + mapInIn[from] = to; + } + + // Traverse the map to find the smallest indices + for(std::map::iterator iter = mapInIn.begin(); + iter != mapInIn.end(); iter++) { + Index end = iter->second; + do { + std::map::const_iterator iter = + mapInIn.find(end); + if (iter != mapInIn.end()) { + end = iter->second; + } else { + break; + } + } while (1); + iter->second = end; + } + + // Create output list + sofa::type::vector mapInOut; + mapInOut.resize(inPt.size()); + + for (Index i=0, newId=0; i::const_iterator iter = + mapInIn.find(i); + if (iter == mapInIn.end()) { + outPt.push_back(inPt[i]); + mapInOut[i] = newId; + imapNode2Node[newId] = i; + newId++; + // Take the position + outMPt[i] = inPt[i]; + outMNorm[i] = inNorm[i]; + } else { + mapInOut[i] = mapInOut[iter->second]; + //imapNode2Node[ mapInOut[iter->second] ] = i; + // Take the position of the point we join with + outMPt[i] = inPt[iter->second]; + outMNorm[i] = inNorm[iter->second]; + } + } + + // Compute new normals by averaging the normals of joined nodes. + if (bHasNormals) { + outNorm.resize(outPt.size()); + for (Index i=0; i(mapInIn, mapInOut, f_input_edges, f_output_edges); + createElements<3>(mapInIn, mapInOut, f_input_triangles, f_output_triangles); + createElements<4>(mapInIn, mapInOut, f_input_quads, f_output_quads); + createElements<4>(mapInIn, mapInOut, f_input_tetrahedra, f_output_tetrahedra); + createElements<8>(mapInIn, mapInOut, f_input_hexahedra, f_output_hexahedra); +} + +template +template +void JoinMeshPoints::createElements( + std::map mapInIn, + sofa::type::vector mapInOut, + const sofa::Data< sofa::type::vector< sofa::type::fixed_array > > &inElements, + sofa::Data< sofa::type::vector< sofa::type::fixed_array > > &outElements) +{ + const sofa::type::vector< sofa::type::fixed_array >& inEle = inElements.getValue(); + + sofa::type::vector< sofa::type::fixed_array >& outEle = *outElements.beginEdit(); + outEle.resize(inEle.size()); + + for (Index i= 0; i& out = outEle[i]; + for (Index j=0; j::const_iterator iter = + mapInIn.find(id); + + bool bJoined = false; + // In -> In mapping (joining nodes) + if (iter != mapInIn.end()) { + id = iter->second; + bJoined = true; + } + + // In -> Out mapping (reindexing nodes) + if (id < f_input_position.getValue().size()) { + out[j] = mapInOut[id]; + } else { + // Invalid node ID, what now? + msg_warning() << "Invalid node ID in elements! " << + id << " is not a valid node ID!" ; + out[j] = 0; + } + + if (N == 3) { + if (bJoined) + imapTriNode2Node[i][out[j]] = inEle[i][j]; + } + + } + } + + outElements.endEdit(); +} + +} // namespace diff --git a/src/Shell/forcefield/TriangularBendingFEMForceField.cpp b/src/Shell/forcefield/TriangularBendingFEMForceField.cpp new file mode 100644 index 0000000..fd1cebf --- /dev/null +++ b/src/Shell/forcefield/TriangularBendingFEMForceField.cpp @@ -0,0 +1,45 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#define SHELL_FORCEFIELD_TRIANGULARBENDINGFEMFORCEFIELD_CPP + +#include +#include +#include +#include + +#include +#include + + +namespace shell::forcefield +{ + +using namespace sofa::defaulttype; + +// Register in the Factory +int TriangularBendingFEMForceFieldClass = sofa::core::RegisterObject("Triangular finite elements with bending") +.add< TriangularBendingFEMForceField >() +; + +template class SOFA_SHELL_API TriangularBendingFEMForceField; + +} // namespace diff --git a/src/Shell/forcefield/TriangularBendingFEMForceField.h b/src/Shell/forcefield/TriangularBendingFEMForceField.h new file mode 100644 index 0000000..8008f8f --- /dev/null +++ b/src/Shell/forcefield/TriangularBendingFEMForceField.h @@ -0,0 +1,295 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include + +#include +#include + +#include +#include + + +namespace shell::forcefield +{ + +using namespace sofa::type; +using sofa::type::vector; +using namespace sofa::core::topology; +using namespace sofa::core::behavior; + +/// This class can be overridden if needed for additionnal storage within template specializations. +template +class TriangularBendingFEMForceFieldInternalData +{ +public: +}; + +template +class TriangularBendingFEMForceField : public sofa::core::behavior::ForceField +{ + public: + SOFA_CLASS(SOFA_TEMPLATE(TriangularBendingFEMForceField,DataTypes), SOFA_TEMPLATE(sofa::core::behavior::ForceField,DataTypes)); + + typedef sofa::core::behavior::ForceField Inherited; + typedef typename DataTypes::VecCoord VecCoord; + typedef typename DataTypes::VecDeriv VecDeriv; + + typedef typename DataTypes::Coord Coord; + typedef typename DataTypes::Deriv Deriv; + typedef typename Coord::value_type Real; + + typedef Vec<3,Real> Vec3; + typedef Vec<2,Real> Vec2; + typedef sofa::type::Quat Quat; + + typedef sofa::Data DataVecCoord; + typedef sofa::Data DataVecDeriv; + + typedef sofa::defaulttype::Vec3Types::VecCoord VecCoordHigh; + + typedef sofa::Index Index; + typedef sofa::core::topology::BaseMeshTopology::Triangle Triangle; + typedef sofa::core::topology::BaseMeshTopology::SeqTriangles SeqTriangles; + + protected: + + typedef Vec<6, Real> Displacement; ///< the displacement vector + typedef Vec<9, Real> DisplacementBending; ///< the displacement vector for bending + typedef Mat<3, 3, Real> MaterialStiffness; ///< the matrix of material stiffness + typedef Mat<6, 3, Real> StrainDisplacement; ///< the strain-displacement matrix + typedef Mat<3, 9, Real> StrainDisplacementBending; + typedef Mat<3, 3, Real > Transformation; ///< matrix for rigid transformations like rotations + typedef Mat<6, 6, Real> StiffnessMatrix; + typedef Mat<9, 9, Real> StiffnessMatrixBending; + typedef Mat<18, 18, Real> StiffnessMatrixGlobalSpace; + + sofa::core::topology::BaseMeshTopology* _topology; + +public: + + class TriangleInformation + { + public: + + // position of B, C vertices in local (in-plane) coordinates (A is naturaly [0,0]) + sofa::type::fixed_array restLocalPositions; + sofa::type::fixed_array restLocalOrientations; + + /// material stiffness matrices of each tetrahedron + MaterialStiffness materialMatrix; + // the strain-displacement matrices vector + StrainDisplacement strainDisplacementMatrix; + // the strain-displacement matrices vector (at Gauss points) + StrainDisplacementBending strainDisplacementMatrix1; + StrainDisplacementBending strainDisplacementMatrix2; + StrainDisplacementBending strainDisplacementMatrix3; + // Indices of each vertex + Index a, b, c; + // Indices in rest shape topology. Normaly is the same as a, b, + // c, but is different if topologyMapper is set. + Index a0, b0, c0; + // Local coordinates + Vec3 localB, localC; + // Transformation rotation; + Quat Qframe; + // Stiffness matrix K = J * M * Jt + StiffnessMatrix stiffnessMatrix; + // Stiffness matrix for bending K = Jt * M * J + StiffnessMatrixBending stiffnessMatrixBending; + + // Surface + Real area; + // Variables needed for drawing the shell + Vec<9, Real> u; // displacement vector + Mat<9, 9, Real> invC; // inverse of C (used in bending mode only) + Vec<9, Real> coefficients; // coefficients Ci computed from tinfo->invC * (tinfo->u + tinfo->u_flat) + Vec <9, Real> u_rest; // difference between the initial position and the flate position to allow the use of an initial deformed shape + + TriangleInformation() { } + + /// Output stream + inline friend std::ostream& operator<< ( std::ostream& os, const TriangleInformation& /*ti*/ ) + { + return os; + } + + /// Input stream + inline friend std::istream& operator>> ( std::istream& in, TriangleInformation& /*ti*/ ) + { + return in; + } + }; + + class TRQSTriangleHandler : public TopologyDataHandler > + { + public: + TRQSTriangleHandler(TriangularBendingFEMForceField* _ff, TriangleData >* _data) : TopologyDataHandler >(_data), ff(_ff) {} + + void applyCreateFunction(unsigned int triangleIndex, TriangleInformation& , + const Triangle & t, + const sofa::type::vector< unsigned int > &, + const sofa::type::vector< double > &); + + protected: + TriangularBendingFEMForceField* ff; + }; + + TriangularBendingFEMForceField(); + + virtual ~TriangularBendingFEMForceField(); + void init() override; + void reinit() override; + void addForce(const sofa::core::MechanicalParams* /*mparams*/, DataVecDeriv& dataF, const DataVecCoord& dataX, const DataVecDeriv& /*dataV*/ ) override ; + void addDForce(const sofa::core::MechanicalParams* /*mparams*/, DataVecDeriv& datadF, const DataVecDeriv& datadX ) override ; + void addKToMatrix(const sofa::core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) override; + void addBToMatrix(sofa::linearalgebra::BaseMatrix * /*mat*/, double /*bFact*/, unsigned int &/*offset*/) override; + double getPotentialEnergy(const VecCoord& x) const; + void handleTopologyChange() override; + + SReal getPotentialEnergy(const sofa::core::MechanicalParams* /*mparams*/, const DataVecCoord& /*x*/) const override { return 0; } + + void draw(const sofa::core::visual::VisualParams* vparams) override; + + sofa::core::topology::BaseMeshTopology* getTopology() {return _topology;} + TriangleData< sofa::type::vector >& getTriangleInfo() {return triangleInfo;} + + sofa::Data d_poisson; + sofa::Data d_young; + sofa::Data d_bending; + sofa::Data d_thickness; + sofa::Data d_membraneRatio; + sofa::Data d_bendingRatio; + sofa::Data d_refineMesh; + sofa::Data d_iterations; + sofa::SingleLink, + sofa::core::topology::BaseMeshTopology, + sofa::BaseLink::FLAG_STOREPATH|sofa::BaseLink::FLAG_STRONGLINK> l_targetTopology; + VecCoordHigh m_targetVertices; + SeqTriangles m_targetTriangles; + + // Allow transition between rest shapes + sofa::SingleLink, + shell::controller::MeshInterpolator, + sofa::BaseLink::FLAG_STOREPATH|sofa::BaseLink::FLAG_STRONGLINK> l_restShape; + + // Indirect rest shape indexing (e.g. for "joining" two meshes) + bool m_mapTopology; + sofa::SingleLink, + shell::engine::JoinMeshPoints, + sofa::BaseLink::FLAG_STOREPATH|sofa::BaseLink::FLAG_STRONGLINK> l_topologyMapper; + + + sofa::core::objectmodel::DataFileName m_exportFilename; + sofa::Data d_exportEveryNbSteps; + sofa::Data d_exportAtBegin; + sofa::Data d_exportAtEnd; + unsigned int m_stepCounter; + + TRQSTriangleHandler* m_triangleHandler; + +protected : + + TriangleData< sofa::type::vector > triangleInfo; + + void computeDisplacement(Displacement &Disp, const VecCoord &x, const Index elementIndex); + void computeDisplacementBending(DisplacementBending &Disp, const VecCoord &x, const Index elementIndex); + void computeStrainDisplacementMatrix(StrainDisplacement &J, const Index elementIndex, const Vec3& b, const Vec3& c); + void computeStrainDisplacementMatrixBending(TriangleInformation *tinfo, const Vec3& b, const Vec3& c); + void tensorFlatPlate(Mat<3, 9, Real>& D, const Vec3 &P); + void computeStiffnessMatrix(StiffnessMatrix &K, const StrainDisplacement &J, const MaterialStiffness &M); + void computeStiffnessMatrixBending(StiffnessMatrixBending &K, TriangleInformation *tinfo); + void computeForce(Displacement &F, const Displacement& D, const Index elementIndex); + void computeForceBending(DisplacementBending &F, const DisplacementBending& D, const Index elementIndex); + + static void TRQSTriangleCreationFunction (unsigned int , void* , TriangleInformation &, const Triangle& , const sofa::type::vector< unsigned int > &, const sofa::type::vector< double >&); + + /// f += Kx where K is the stiffness matrix and x a displacement + virtual void applyStiffness(VecDeriv& f, const VecDeriv& dx, const Index elementIndex, const double kFactor); + virtual void computeMaterialStiffness(const int i); + + void initTriangleOnce(const int i, const Index&a, const Index&b, const Index&c); + void initTriangle(const int i); + void computeRotation(Quat &Qframe, const VecCoord &p, const Index &a, const Index &b, const Index &c); + void accumulateForce(VecDeriv& f, const VecCoord & p, const Index elementIndex); + + void convertStiffnessMatrixToGlobalSpace(StiffnessMatrixGlobalSpace &K_gs, TriangleInformation *tinfo); + + void refineCoarseMeshToTarget(void); + void subdivide(const Vec3& a, const Vec3& b, const Vec3& c, sofa::type::vector &subVertices, SeqTriangles &subTriangles); + void addVertexAndFindIndex(sofa::type::vector &subVertices, const Vec3 &vertex, int &index); + void movePoint(Vec3& pointToMove); + void findClosestGravityPoints(const Vec3& point, sofa::type::vector& listClosestPoints); + + void handleEvent(sofa::core::objectmodel::Event *event) override; + + + Quat qDiff(Quat a, const Quat& b) + { + // If the axes are not oriented in the same direction, flip the axis and angle of a to get the same convention than b + if (a[0]*b[0]+a[1]*b[1]+a[2]*b[2]+a[3]*b[3]<0) + { + a[0] = -a[0]; + a[1] = -a[1]; + a[2] = -a[2]; + a[3] = -a[3]; + } + Quat q = b.inverse() * a; + return q; + } + + Quat qDiffZ(const Quat& vertex, const Quat& Qframe) + { + // dQ is the quaternion that embodies the rotation between the z axis of the vertex and the z axis of the local triangle's frame (in local space) + Quat dQ; + + // u = z axis of the triangle's frame + Vec3d u(0,0,1); + + // v = z axis of the vertex's frame is expressed into world space + Vec3d v = vertex.rotate(Vec3d(0.0, 0.0, 1.0)); + // v0 = v expressed into local triangle's frame + Vec3d v0 = Qframe.rotate(v); + + // Axis of rotation between the 2 vectors u and v lies into the plan of the 2 vectors + Vec3d axis = cross(u, v0); + // Shortest angle between the 2 vectors + double angle = acos(dot(u, v0)); + + // Quaternion associated to this axis and this angle + if (fabs(angle)>1e-6) + { + dQ.axisToQuat(axis,angle); + } + else + { + dQ = Quat(0,0,0,1); + } + + return dQ; + } +}; + +} // namespace diff --git a/src/Shell/forcefield/TriangularBendingFEMForceField.inl b/src/Shell/forcefield/TriangularBendingFEMForceField.inl new file mode 100644 index 0000000..122957d --- /dev/null +++ b/src/Shell/forcefield/TriangularBendingFEMForceField.inl @@ -0,0 +1,1262 @@ +/****************************************************************************** +* SOFA, Simulation Open-Framework Architecture * +* (c) 2006 INRIA, USTL, UJF, CNRS, MGH * +* * +* This program is free software; you can redistribute it and/or modify it * +* under the terms of the GNU Lesser General Public License as published by * +* the Free Software Foundation; either version 2.1 of the License, or (at * +* your option) any later version. * +* * +* This program is distributed in the hope that it will be useful, but WITHOUT * +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * +* FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License * +* for more details. * +* * +* You should have received a copy of the GNU Lesser General Public License * +* along with this program. If not, see . * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#pragma once + +#include +#include +#include // for reading the file +#include //for debugging +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +#ifdef _WIN32 +#include +#endif + +namespace shell::forcefield +{ + using namespace sofa::type; + // using namespace sofa::component::topology; + + + +// -------------------------------------------------------------------------------------- +// --- Topology Creation/Destruction functions +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::TRQSTriangleHandler::applyCreateFunction(unsigned int triangleIndex, TriangleInformation &, const Triangle &t, const sofa::type::vector &, const sofa::type::vector &) +{ + if (ff) + { + ff->initTriangleOnce(triangleIndex, t[0], t[1], t[2]); + ff->initTriangle(triangleIndex); + ff->computeMaterialStiffness(triangleIndex); + } +} + + +// -------------------------------------------------------------------------------------- +// --- +// -------------------------------------------------------------------------------------- +template +TriangularBendingFEMForceField::TriangularBendingFEMForceField() +: d_poisson(initData(&d_poisson,(Real)0.45,"poissonRatio","Poisson ratio in Hooke's law")) +, d_young(initData(&d_young,(Real)3000.,"youngModulus","Young modulus in Hooke's law")) +, d_bending(initData(&d_bending,false,"bending","Adds bending")) +, d_thickness(initData(&d_thickness,(Real)0.1,"thickness","Thickness of the plates")) +, d_membraneRatio(initData(&d_membraneRatio,(Real)1.0,"membraneRatio","In plane forces ratio")) +, d_bendingRatio(initData(&d_bendingRatio,(Real)1.0,"bendingRatio","Bending forces ratio")) +, d_refineMesh(initData(&d_refineMesh, false, "refineMesh","Hierarchical refinement of the mesh")) +, d_iterations(initData(&d_iterations,(int)0,"iterations","Iterations for refinement")) +, l_targetTopology(initLink("targetTopology","Targeted high resolution topology")) +, l_restShape(initLink("restShape","MeshInterpolator component for variable rest shape")) +, m_mapTopology(false) +, l_topologyMapper(initLink("topologyMapper","Component supplying different topology for the rest shape")) +, m_exportFilename(initData(&m_exportFilename, "exportFilename", "file name to export coefficients into")) +, d_exportEveryNbSteps(initData(&d_exportEveryNbSteps, (unsigned int)0, "exportEveryNumberOfSteps", "export file only at specified number of steps (0=disable)")) +, d_exportAtBegin(initData(&d_exportAtBegin, false, "exportAtBegin", "export file at the initialization")) +, d_exportAtEnd(initData(&d_exportAtEnd, false, "exportAtEnd", "export file when the simulation is finished")) +, m_stepCounter(0) +, triangleInfo(initData(&triangleInfo, "triangleInfo", "Internal triangle data")) +{ + m_triangleHandler = new TRQSTriangleHandler(this, &triangleInfo); +} + +// -------------------------------------------------------------------------------------- +// --- +// -------------------------------------------------------------------------------------- +template void TriangularBendingFEMForceField::handleTopologyChange() +{ +} + +// -------------------------------------------------------------------------------------- +// --- +// -------------------------------------------------------------------------------------- +template +TriangularBendingFEMForceField::~TriangularBendingFEMForceField() +{ + if(m_triangleHandler) delete m_triangleHandler; +} + +// -------------------------------------------------------------------------------------- +// --- Initialization stage +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::init() +{ + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Valid); + this->Inherited::init(); + + _topology = this->getContext()->getMeshTopology(); + + if (_topology->getNbTriangles()==0) + { + msg_error() << "TriangularBendingFEMForceField: object must have a Triangular Set Topology."; + this->d_componentState.setValue(sofa::core::objectmodel::ComponentState::Invalid); + return; + } + + // Create specific handler for TriangleData + triangleInfo.createTopologyHandler(_topology); + + reinit(); + + if (d_refineMesh.getValue()) + { + sofa::core::topology::BaseMeshTopology* _topologyTarget = l_targetTopology.get(); + + if (_topologyTarget) + { + MechanicalState* mStateTarget = dynamic_cast*> (_topologyTarget->getContext()->getMechanicalState()); + if (mStateTarget) + { + m_targetTriangles = _topologyTarget->getTriangles(); + m_targetVertices = mStateTarget->read(sofa::core::vec_id::read_access::position)->getValue(); + } + else + { + msg_warning() << "No mechanical state for target high resolution topology" ; + return; + } + } + else + { + msg_warning() << "No target high resolution mesh found"; + return; + } + + // Run procedure for shell remeshing + refineCoarseMeshToTarget(); + } +} + + +template +void TriangularBendingFEMForceField::refineCoarseMeshToTarget(void) +{ + msg_info() << "Refining a mesh of " << _topology->getNbTriangles() << " triangles towards a target surface of " << m_targetTriangles.size() << " triangles."; + + // List of vertices + const VecCoord& x = this->mstate->read(sofa::core::vec_id::read_access::position)->getValue(); + // List of triangles + const SeqTriangles triangles = _topology->getTriangles(); + + // Creates new mesh + sofa::type::vector subVertices; + SeqTriangles subTriangles; + + // Initialises list of subvertices and triangles + for (unsigned int i=0; i +void TriangularBendingFEMForceField::subdivide(const Vec3& a, const Vec3& b, const Vec3& c, sofa::type::vector &subVertices, SeqTriangles &subTriangles) +{ + // Global coordinates + Vec3 mAB, mAC, mBC; + mAB = (a+b)/2; + mAC = (a+c)/2; + mBC = (b+c)/2; + + // Adds vertex if we deal with a new point + int indexAB, indexAC, indexBC; + addVertexAndFindIndex(subVertices, mAB, indexAB); + addVertexAndFindIndex(subVertices, mAC, indexAC); + addVertexAndFindIndex(subVertices, mBC, indexBC); + + // Finds index of the 3 original vertices + int indexA, indexB, indexC; + addVertexAndFindIndex(subVertices, a, indexA); + addVertexAndFindIndex(subVertices, b, indexB); + addVertexAndFindIndex(subVertices, c, indexC); + + // Adds the 4 subdivided triangles to the list + subTriangles.push_back(Triangle(indexA, indexAB, indexAC)); + subTriangles.push_back(Triangle(indexAB, indexB, indexBC)); + subTriangles.push_back(Triangle(indexAC, indexBC, indexC)); + subTriangles.push_back(Triangle(indexBC, indexAC, indexAB)); +} + + +// -------------------------------------------------------------------------------------- +// Adds a vertex if it is not already in the list +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::addVertexAndFindIndex(sofa::type::vector &subVertices, const Vec3 &vertex, int &index) +{ + bool alreadyHere = false; + + for (unsigned int v=0; v +void TriangularBendingFEMForceField::movePoint(Vec3& pointToMove) +{ + sofa::type::vector listClosestPoints; + findClosestGravityPoints(pointToMove, listClosestPoints); + pointToMove = (listClosestPoints[0]+listClosestPoints[1]+listClosestPoints[2])/3; +} + + +// -------------------------------------------------------------------------------------- +// Finds the list of the 3 closest gravity points of targeted surface +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::findClosestGravityPoints(const Vec3& point, sofa::type::vector& listClosestPoints) +{ + std::multimap closestTrianglesData; + + for (unsigned int t=0; t(1.0,Vec3()); + closestTrianglesData.insert( std::make_pair(distance,G)); + } + + // Returns the 3 closest points + int count = 0; + typename std::multimap::iterator it; + for (it = closestTrianglesData.begin(); it!=closestTrianglesData.end(); it++) + { + if (count < 3) + { + listClosestPoints.push_back((*it).second); + } + count++; + } + +} + + +// -------------------------------------------------------------------------------------- +// --- Re-initialization (called when we change a parameter through the GUI) +// -------------------------------------------------------------------------------------- +template void TriangularBendingFEMForceField::reinit() +{ + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + + if (l_topologyMapper.get() != nullptr) { + + shell::engine::JoinMeshPoints* jmp = l_topologyMapper.get(); + if (jmp->f_output_triangles.getValue().size() == 0) + { + msg_warning() << "Mapped topology must be triangular. No triangles found." ; + } else { + m_mapTopology = true; + } + } + + if (l_restShape.get() != nullptr) { + // Listen for MeshChangedEvent + *this->f_listening.beginEdit() = true; + this->f_listening.endEdit(); + + // Check if there is same number of nodes + if (!m_mapTopology) { + if (l_restShape.get()->f_position.getValue().size() != + this->mstate->read(sofa::core::vec_id::read_access::position)->getValue().size()) { + msg_warning() << "Different number of nodes in rest shape and mechanical state." ; + } + } else if (l_restShape.get()->f_position.getValue().size() != + l_topologyMapper.get()->f_input_position.getValue().size()) { + msg_warning() << "Different number of nodes in rest shape and (original) mapped topology." ; + } + } + + /// Prepare to store info in the triangle array + triangleInf.resize(_topology->getNbTriangles()); + + for (sofa::Index i=0; i<_topology->getNbTriangles(); ++i) + { + m_triangleHandler->applyCreateFunction(i, triangleInf[i], _topology->getTriangle(i), (const sofa::type::vector< unsigned int > )0, (const sofa::type::vector< double >)0); + } + + triangleInfo.endEdit(); +} + + +// -------------------------------------------------------------------------------------- +// --- +// -------------------------------------------------------------------------------------- +template +double TriangularBendingFEMForceField::getPotentialEnergy(const VecCoord& /*x*/) const +{ + msg_warning() << "TriangularBendingFEMForceField::getPotentialEnergy is not implemented."; + return 0; +} + + +// -------------------------------------------------------------------------------------- +// Computes the quaternion that embodies the rotation from triangle to world +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::computeRotation(Quat& Qframe, const VecCoord &x, const Index &a, const Index &b, const Index &c) +{ + // First vector on first edge + // Second vector in the plane of the two first edges + // Third vector orthogonal to first and second + + Vec3 edgex = x[b].getCenter() - x[a].getCenter(); + edgex.normalize(); + + Vec3 edgey = x[c].getCenter() - x[a].getCenter(); + edgey.normalize(); + + Vec3 edgez; + edgez = cross(edgex, edgey); + edgez.normalize(); + + edgey = cross(edgez, edgex); + edgey.normalize(); + + Transformation R; + R[0][0] = edgex[0]; + R[0][1] = edgex[1]; + R[0][2] = edgex[2]; + R[1][0] = edgey[0]; + R[1][1] = edgey[1]; + R[1][2] = edgey[2]; + R[2][0] = edgez[0]; + R[2][1] = edgez[1]; + R[2][2] = edgez[2]; + + Qframe.fromMatrix(R); +} + + +// -------------------------------------------------------------------------------------- +// --- Initialization of a triangle, is called *only* once +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::initTriangleOnce(const int i, const Index&a, const Index&b, const Index&c) +{ + + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + TriangleInformation *tinfo = &triangleInf[i]; + + // Store indices of each vertex + tinfo->a = a; + tinfo->b = b; + tinfo->c = c; + + Index a0=a, b0=b, c0=c; + if (m_mapTopology) { + shell::engine::JoinMeshPoints* jmp = l_topologyMapper.get(); + + // Get indices in original topology + a0 = jmp->getSrcNodeFromTri(i, a0); + b0 = jmp->getSrcNodeFromTri(i, b0); + c0 = jmp->getSrcNodeFromTri(i, c0); + } + + tinfo->a0 = a0; + tinfo->b0 = b0; + tinfo->c0 = c0; + + triangleInfo.endEdit(); +} + +// -------------------------------------------------------------------------------------- +// --- Initialization of a triangle. Can be called more than once if there is a +// --- changing rest shape. +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::initTriangle(const int i) +{ + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + TriangleInformation *tinfo = &triangleInf[i]; + + const Index a = tinfo->a; + const Index b = tinfo->b; + const Index c = tinfo->c; + + const Index a0 = tinfo->a0; + const Index b0 = tinfo->b0; + const Index c0 = tinfo->c0; + + // Gets vertices of rest and initial positions respectively + const VecCoord& x0 = (l_restShape.get() != nullptr) + // if having changing rest shape take it + ? l_restShape.get()->f_position.getValue() + : (m_mapTopology + // if rest shape is fixed but we have mapped topology use it + ? l_topologyMapper.get()->f_input_position.getValue() + // otherwise just take rest shape in mechanical state + : this->mstate->read(sofa::core::vec_id::read_access::position)->getValue() + ); + const VecCoord& x = this->mstate->read(sofa::core::vec_id::read_access::position)->getValue(); + + // Rotation from triangle to world at rest and initial positions (respectively) + Quat Qframe0, Qframe; + computeRotation(Qframe0, x0, a0, b0, c0 ); + computeRotation(Qframe, x, a, b, c ); + tinfo->Qframe = Qframe; + + // The positions of each point is expressed into the local frame at rest position + tinfo->restLocalPositions[0] = Qframe0.rotate(x0[b0].getCenter() - x0[a0].getCenter()); + tinfo->restLocalPositions[1] = Qframe0.rotate(x0[c0].getCenter() - x0[a0].getCenter()); + + if (d_bending.getValue()) + { + // Computes inverse of C for initial position (in case of the latter is different than the rest_position) + tinfo->localB = Qframe.rotate(x[b].getCenter()-x[a].getCenter()); + tinfo->localC = Qframe.rotate(x[c].getCenter()-x[a].getCenter()); + computeStrainDisplacementMatrixBending(tinfo, tinfo->localB, tinfo->localC); + + // Computes triangles' surface + StrainDisplacement J; + computeStrainDisplacementMatrix(J, i, tinfo->localB, tinfo->localC); + + // Local rest orientations (Evaluates the difference between the rest position and the flat position to allow the use of a deformed rest shape) + tinfo->restLocalOrientations[0] = qDiffZ(x0[a0].getOrientation(), Qframe0); + tinfo->restLocalOrientations[1] = qDiffZ(x0[b0].getOrientation(), Qframe0); + tinfo->restLocalOrientations[2] = qDiffZ(x0[c0].getOrientation(), Qframe0); + + // Creates a vector u_rest matching the difference between rest and flat positions + tinfo->u_rest.clear(); + tinfo->u_rest[1] = tinfo->restLocalOrientations[0].toEulerVector()[0]; + tinfo->u_rest[2] = tinfo->restLocalOrientations[0].toEulerVector()[1]; + tinfo->u_rest[4] = tinfo->restLocalOrientations[1].toEulerVector()[0]; + tinfo->u_rest[5] = tinfo->restLocalOrientations[1].toEulerVector()[1]; + tinfo->u_rest[7] = tinfo->restLocalOrientations[2].toEulerVector()[0]; + tinfo->u_rest[8] = tinfo->restLocalOrientations[2].toEulerVector()[1]; + + // Computes vector displacement between initial position and rest positions (actual displacements that define the amount of stress within the structure) + DisplacementBending Disp_bending; + computeDisplacementBending(Disp_bending, x, i); + + } + + triangleInfo.endEdit(); +} + +// -------------------------------------------------------------------------------------- +// --- +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::applyStiffness(VecDeriv& v, const VecDeriv& dx, const Index elementIndex, const double kFactor) +{ + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + TriangleInformation *tinfo = &triangleInf[elementIndex]; + + // Get the indices of the 3 vertices for the current triangle + const Index& a = tinfo->a; + const Index& b = tinfo->b; + const Index& c = tinfo->c; + + // Computes displacements + Displacement Disp; + Vec3 x_a, x_b, x_c; + x_a = tinfo->Qframe.rotate(getVCenter(dx[a])); + Disp[0] = x_a[0]; + Disp[1] = x_a[1]; + + x_b = tinfo->Qframe.rotate(getVCenter(dx[b])); + Disp[2] = x_b[0]; + Disp[3] = x_b[1]; + + x_c = tinfo->Qframe.rotate(getVCenter(dx[c])); + Disp[4] = x_c[0]; + Disp[5] = x_c[1]; + + // Compute dF + Displacement dF; + dF = tinfo->stiffnessMatrix * Disp; + + // Transfer into global frame + getVCenter(v[a]) += tinfo->Qframe.inverseRotate(Vec3(-dF[0], -dF[1], 0)) * kFactor; + getVCenter(v[b]) += tinfo->Qframe.inverseRotate(Vec3(-dF[2], -dF[3], 0)) * kFactor; + getVCenter(v[c]) += tinfo->Qframe.inverseRotate(Vec3(-dF[4], -dF[5], 0)) * kFactor; + + // If bending is requested + if (d_bending.getValue()) + { + // Bending displacements + DisplacementBending Disp_bending; + Vec3 u; + u = tinfo->Qframe.rotate(getVOrientation(dx[a])); + Disp_bending[0] = x_a[2]; + Disp_bending[1] = u[0]; + Disp_bending[2] = u[1]; + + u = tinfo->Qframe.rotate(getVOrientation(dx[b])); + Disp_bending[3] = x_b[2]; + Disp_bending[4] = u[0]; + Disp_bending[5] = u[1]; + + u = tinfo->Qframe.rotate(getVOrientation(dx[c])); + Disp_bending[6] = x_c[2]; + Disp_bending[7] = u[0]; + Disp_bending[8] = u[1]; + + // Compute dF + DisplacementBending dF_bending; + dF_bending = tinfo->stiffnessMatrixBending * Disp_bending; + + // Go back into global frame + Vec3 fa1, fa2, fb1, fb2, fc1, fc2; + fa1 = tinfo->Qframe.inverseRotate(Vec3(0.0, 0.0, dF_bending[0])); + fa2 = tinfo->Qframe.inverseRotate(Vec3(dF_bending[1], dF_bending[2], 0.0)); + + fb1 = tinfo->Qframe.inverseRotate(Vec3(0.0, 0.0, dF_bending[3])); + fb2 = tinfo->Qframe.inverseRotate(Vec3(dF_bending[4], dF_bending[5], 0.0)); + + fc1 = tinfo->Qframe.inverseRotate(Vec3(0.0, 0.0, dF_bending[6])); + fc2 = tinfo->Qframe.inverseRotate(Vec3(dF_bending[7], dF_bending[8], 0.0)); + + v[a] += Deriv(-fa1, -fa2) * kFactor; + v[b] += Deriv(-fb1, -fb2) * kFactor; + v[c] += Deriv(-fc1, -fc2) * kFactor; + } + + + triangleInfo.endEdit(); +} + +// ------------------------------------------------------------------------------------------------------------- +// --- Compute displacement vector D as the difference between current current position 'p' and initial position +// --- expressed in the co-rotational frame of reference +// ------------------------------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::computeDisplacement(Displacement &Disp, const VecCoord &x, const Index elementIndex) +{ + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + TriangleInformation *tinfo = &triangleInf[elementIndex]; + + Index a = tinfo->a; + Index b = tinfo->b; + Index c = tinfo->c; + + // Compute local postions of vertices b and c in the co-rotational frame (a is always (0,0,0)) + tinfo->localB = tinfo->Qframe.rotate(x[b].getCenter()-x[a].getCenter()); + tinfo->localC = tinfo->Qframe.rotate(x[c].getCenter()-x[a].getCenter()); + + // In-plane local displacements + Vec3 uAB, uAC; + uAB = tinfo->localB - tinfo->restLocalPositions[0]; + uAC = tinfo->localC - tinfo->restLocalPositions[1]; + + Disp[0] = 0; + Disp[1] = 0; + Disp[2] = uAB[0]; + Disp[3] = 0; + Disp[4] = uAC[0]; + Disp[5] = uAC[1]; + + triangleInfo.endEdit(); +} + + +// ------------------------------------------------------------------------------------------------------------- +// --- Compute bending displacement vector D as the difference between current current position 'p' and initial position +// --- expressed in the co-rotational frame of reference +// ------------------------------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::computeDisplacementBending(DisplacementBending &Disp, const VecCoord &x, const Index elementIndex) +{ + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + TriangleInformation *tinfo = &triangleInf[elementIndex]; + + Index a = tinfo->a; + Index b = tinfo->b; + Index c = tinfo->c; + + Quat dQA = qDiffZ(x[a].getOrientation(), tinfo->Qframe); // Rotation of axis Z from Qframe to axis Z from x[a].getOrientation() + Quat dQB = qDiffZ(x[b].getOrientation(), tinfo->Qframe); // Rotation of axis Z from Qframe to axis Z from x[b].getOrientation() + Quat dQC = qDiffZ(x[c].getOrientation(), tinfo->Qframe); // Rotation of axis Z from Qframe to axis Z from x[c].getOrientation() + + // Difference between the current and the rest orientations yields displacement (into the triangle's frame) + dQA = qDiff(tinfo->restLocalOrientations[0].inverse(), dQA.inverse()); + dQB = qDiff(tinfo->restLocalOrientations[1].inverse(), dQB.inverse()); + dQC = qDiff(tinfo->restLocalOrientations[2].inverse(), dQC.inverse()); + + // Takes the Euler vector to get the rotation's axis + Vec3 rA, rB, rC; + rA = dQA.toEulerVector(); + rB = dQB.toEulerVector(); + rC = dQC.toEulerVector(); + + // Writes the computed displacements + Disp[0] = 0; // z displacement in A + Disp[1] = rA[0]; // x rotation in A + Disp[2] = rA[1]; // y rotation in A + + Disp[3] = 0; // z displacement in B + Disp[4] = rB[0]; // x rotation in B + Disp[5] = rB[1]; // y rotation in B + + Disp[6] = 0; // z displacement in C + Disp[7] = rC[0]; // x rotation in C + Disp[8] = rC[1]; // y rotation in C + + // Stores the vector u of displacements (used by the mechanical mapping for rendering) + tinfo->u = Disp; + + triangleInfo.endEdit(); +} + +// ------------------------------------------------------------------------------------------------------------ +// --- Compute the strain-displacement matrix where (a, b, c) are the local coordinates of the 3 nodes of a triangle +// ------------------------------------------------------------------------------------------------------------ +template +void TriangularBendingFEMForceField::computeStrainDisplacementMatrix(StrainDisplacement &J, const Index elementIndex, const Vec3& b, const Vec3& c) +{ + Real determinant; + determinant = b[0] * c[1]; + + Real x13 = -c[0] / determinant; // since a=(0,0) + Real x21 = b[0] / determinant; // since a=(0,0) + Real x32 = (c[0]-b[0]) / determinant; + Real y12 = 0; // since a=(0,0) and b[1] = 0 + Real y23 = -c[1] / determinant; // since a=(0,0) and b[1] = 0 + Real y31 = c[1] / determinant; // since a=(0,0) + + J[0][0] = y23; + J[0][1] = 0; + J[0][2] = x32; + + J[1][0] = 0; + J[1][1] = x32; + J[1][2] = y23; + + J[2][0] = y31; + J[2][1] = 0; + J[2][2] = x13; + + J[3][0] = 0; + J[3][1] = x13; + J[3][2] = y31; + + J[4][0] = y12; + J[4][1] = 0; + J[4][2] = x21; + + J[5][0] = 0; + J[5][1] = x21; + J[5][2] = y12; + + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + triangleInf[elementIndex].area = 0.5*determinant; + triangleInfo.endEdit(); + +} + + +// ------------------------------------------------------------------------------------------------------------ +// --- Compute the bending strain-displacement matrix where (a, b, c) are the coordinates of the 3 nodes of a triangle +// ------------------------------------------------------------------------------------------------------------ +template +void TriangularBendingFEMForceField::computeStrainDisplacementMatrixBending(TriangleInformation *tinfo, const Vec3& b, const Vec3& c) +{ + // Computation of the inverse of matrix C. Its inverse gives the coefficients c1, c2, ..., c9 of the deflection function given by: + // Uz = c1 + c2*x+ c3*y + c4*x^2 + c5*x*y + c6*y^2 + c7*x^3 + c8*(x*y^2 + x^2*y) + c9*y^3 + // Source: Tocher's deflection function presented by Przemieniecki + // Corrected deflection to get a symmetrical deformation: Uz = c1 + c2*x+ c3*y + c4*x^2 + c5*x*y + c6*y^2 + c7*x^3 + c8*x*y^2 + c9*y^3 + + Mat<9, 9, Real> C; + C.clear(); + + // Corrected + C(0,0) = 1; + C(1,2) = 1; + C(2,1) = -1; + C(3,0) = 1; C(3,1) = b[0]; C(3,3) = b[0]*b[0]; C(3,6) = b[0]*b[0]*b[0]; + C(4,2) = 1; C(4,4) = b[0]; + C(5,1) = -1; C(5,3) = -2*b[0]; C(5,6) = -3*b[0]*b[0]; + C(6,0) = 1; C(6,1) = c[0]; C(6,2) = c[1]; C(6,3) = c[0]*c[0]; C(6,4) = c[0]*c[1]; C(6,5) = c[1]*c[1]; + C(6,6) = c[0]*c[0]*c[0]; C(6,7) = c[0]*c[1]*c[1]; C(6,8) = c[1]*c[1]*c[1]; + C(7,2) = 1; C(7,4) = c[0]; C(7,5) = 2*c[1]; C(7,7) = 2*c[0]*c[1]; C(7,8) = 3*c[1]*c[1]; + C(8,1) = -1; C(8,3) = -2*c[0]; C(8,4) = -c[1]; C(8,6) = -3*c[0]*c[0]; C(8,7) = -c[1]*c[1]; + + // Inverse of C + Mat<9, 9, Real> invC; + invC.invert(C); + tinfo->invC = invC; + + // Calculation of the 3 Gauss points taken in the centre of each edge + Vec3 gaussPoint1 = b*0.5; // a is (0, 0, 0) + Vec3 gaussPoint2 = (b+c)*0.5; + Vec3 gaussPoint3 = c*0.5; // a is (0, 0, 0) + + // Retrieves the strain tensor used in flat-plate theory at each Gauss point + Mat<3, 9, Real> D1, D2, D3; + tensorFlatPlate(D1, gaussPoint1); + tensorFlatPlate(D2, gaussPoint2); + tensorFlatPlate(D3, gaussPoint3); + + // Compute strain-displacement matrix + tinfo->strainDisplacementMatrix1 = D1 * invC; + tinfo->strainDisplacementMatrix2 = D2 * invC; + tinfo->strainDisplacementMatrix3 = D3 * invC; +} + + +// -------------------------------------------------------------------------------------------------------- +// --- Computes the strain tensor used in flat-plate theory in a given point +// -------------------------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::tensorFlatPlate(Mat<3, 9, Real>& D, const Vec3 &P) +{ + // Flat-plat theory gives: + // e = D * c with + // [ 0 0 0 2 0 0 6x 2y 0 ] + // D = -z | 0 0 0 0 0 2 0 2x 6y | + // [ 0 0 0 0 2 0 0 4(x+y) 0 ] + // where e is the strain vector and c the coefficient vector of the deflection function + // + // CORRECTED: + // [ 0 0 0 2 0 0 6x 0 0 ] + // D = -z | 0 0 0 0 0 2 0 2x 6y | + // [ 0 0 0 0 2 0 0 4y 0 ] + + // Corrected + D.clear(); + D(0,3) = 2; D(0,6) = 6*P[0]; + D(1,5) = 2; D(1,7) = 2*P[0]; D(1,8) = 6*P[1]; + D(2,4) = 2; D(2,7) = 4*P[1]; +} + + +// ---------------------------------------------------------------------------------------------------------------------- +// --- Compute the stiffness matrix K = J * M * Jt where J is the strain-displacement matrix and M the material matrix +// ---------------------------------------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::computeStiffnessMatrix(StiffnessMatrix &K, const StrainDisplacement &J, const MaterialStiffness &M) +{ + Mat<3,6,Real> Jt; + Jt.transpose(J); + + K = J * M * Jt; +} + + +// ---------------------------------------------------------------------------------------------------------------------- +// --- Compute the stiffness matrix for bending K = J * M * Jt where J is the strain-displacement matrix and M the material matrix +// ---------------------------------------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::computeStiffnessMatrixBending(StiffnessMatrixBending &K, TriangleInformation *tinfo) +{ + Mat<9, 3, Real> J1t, J2t, J3t; + J1t.transpose(tinfo->strainDisplacementMatrix1); + J2t.transpose(tinfo->strainDisplacementMatrix2); + J3t.transpose(tinfo->strainDisplacementMatrix3); + + // K = (surface/3)*(t^3)*(J1t*material*J1 + J2t*material*J2 + J3t*material*J3) + K = J1t * tinfo->materialMatrix * tinfo->strainDisplacementMatrix1 + + J2t * tinfo->materialMatrix * tinfo->strainDisplacementMatrix2 + + J3t * tinfo->materialMatrix * tinfo->strainDisplacementMatrix3; +} + +// -------------------------------------------------------------------------------------- +// --- Compute material stiffness +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::computeMaterialStiffness(const int i) +{ + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + + TriangleInformation *tinfo = &triangleInf[i]; + + tinfo->materialMatrix[0][0] = 1; + tinfo->materialMatrix[0][1] = d_poisson.getValue(); + tinfo->materialMatrix[0][2] = 0; + tinfo->materialMatrix[1][0] = d_poisson.getValue(); + tinfo->materialMatrix[1][1] = 1; + tinfo->materialMatrix[1][2] = 0; + tinfo->materialMatrix[2][0] = 0; + tinfo->materialMatrix[2][1] = 0; + tinfo->materialMatrix[2][2] = 0.5f * (1 - d_poisson.getValue()); + + tinfo->materialMatrix *= (d_young.getValue() / (12 * (1 - d_poisson.getValue() * d_poisson.getValue()))); + + triangleInfo.endEdit(); +} + + +// -------------------------------------------------------------------------------------- +// --- Compute force F = J * material * Jt * u +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::computeForce(Displacement &F, const Displacement& D, const Index elementIndex) +{ + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + TriangleInformation *tinfo = &triangleInf[elementIndex]; + + // Compute strain-displacement matrix J + StrainDisplacement J; + computeStrainDisplacementMatrix(J, elementIndex, tinfo->localB, tinfo->localC); + tinfo->strainDisplacementMatrix = J; + + // Compute stiffness matrix K = J*material*Jt + StiffnessMatrix K; + computeStiffnessMatrix(K, J, tinfo->materialMatrix); + tinfo->stiffnessMatrix = K; + + // Compute forces + F = K * D; + + triangleInfo.endEdit(); +} + + +// -------------------------------------------------------------------------------------- +// --- Compute force F = Jt * material * J * u +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::computeForceBending(DisplacementBending &F_bending, const DisplacementBending& D_bending, const Index elementIndex) +{ + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + TriangleInformation *tinfo = &triangleInf[elementIndex]; + + // Compute strain-displacement matrix J + computeStrainDisplacementMatrixBending(tinfo, tinfo->localB, tinfo->localC); + + // Compute stiffness matrix K = Jt * material * J + StiffnessMatrixBending K_bending; + computeStiffnessMatrixBending(K_bending, tinfo); + Real t = d_thickness.getValue(); + K_bending *= t*t*t*(triangleInf[elementIndex].area)/3; + tinfo->stiffnessMatrixBending = K_bending; + + // Compute forces + F_bending = K_bending * D_bending; + + triangleInfo.endEdit(); +} + +// -------------------------------------------------------------------------------------- +// --- +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::accumulateForce(VecDeriv &f, const VecCoord &x, const Index elementIndex) +{ + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + TriangleInformation *tinfo = &triangleInf[elementIndex]; + + // Get the indices of the 3 vertices for the current triangle + const Index& a = tinfo->a; + const Index& b = tinfo->b; + const Index& c = tinfo->c; + + // Compute the quaternion that embodies the rotation between the triangle and world frames (co-rotational method) + Quat Qframe; + computeRotation(Qframe, x, a, b, c); + tinfo->Qframe = Qframe; + + // Compute in-plane displacement into the triangle's frame + Displacement D; + computeDisplacement(D, x, elementIndex); + + // Compute in-plane forces on this element (in the co-rotational space) + Displacement F; + computeForce(F, D, elementIndex); + + // Transform forces back into global reference frame + getVCenter(f[a]) -= tinfo->Qframe.inverseRotate(Vec3(F[0], F[1], 0)); + getVCenter(f[b]) -= tinfo->Qframe.inverseRotate(Vec3(F[2], F[3], 0)); + getVCenter(f[c]) -= tinfo->Qframe.inverseRotate(Vec3(F[4], F[5], 0)); + + if (d_bending.getValue()) + { + // Compute bending displacement for bending into the triangle's frame + DisplacementBending D_bending; + computeDisplacementBending(D_bending, x, elementIndex); + + // Compute bending forces on this element (in the co-rotational space) + DisplacementBending F_bending; + computeForceBending(F_bending, D_bending, elementIndex); + + // Transform forces back into global reference frame + Vec3 fa1 = tinfo->Qframe.inverseRotate(Vec3(0.0, 0.0, F_bending[0])); + Vec3 fa2 = tinfo->Qframe.inverseRotate(Vec3(F_bending[1], F_bending[2], 0.0)); + + Vec3 fb1 = tinfo->Qframe.inverseRotate(Vec3(0.0, 0.0, F_bending[3])); + Vec3 fb2 = tinfo->Qframe.inverseRotate(Vec3(F_bending[4], F_bending[5], 0.0)); + + Vec3 fc1 = tinfo->Qframe.inverseRotate(Vec3(0.0, 0.0, F_bending[6])); + Vec3 fc2 = tinfo->Qframe.inverseRotate(Vec3(F_bending[7], F_bending[8], 0.0)); + + f[a] += Deriv(-fa1, -fa2); + f[b] += Deriv(-fb1, -fb2); + f[c] += Deriv(-fc1, -fc2); + } + + triangleInfo.endEdit(); +} + + +// -------------------------------------------------------------------------------------- +// --- +// -------------------------------------------------------------------------------------- +template +void TriangularBendingFEMForceField::addForce(const sofa::core::MechanicalParams* /*mparams*/, DataVecDeriv& dataF, const DataVecCoord& dataX, const DataVecDeriv& /*dataV*/ ) +{ + VecDeriv& f = *(dataF.beginEdit()); + const VecCoord& p = dataX.getValue() ; + + int nbTriangles=_topology->getNbTriangles(); + f.resize(p.size()); + + for (int i=0; i +void TriangularBendingFEMForceField::addDForce(const sofa::core::MechanicalParams* mparams, DataVecDeriv& datadF, const DataVecDeriv& datadX ) +{ + VecDeriv& df = *(datadF.beginEdit()); + const VecDeriv& dp = datadX.getValue() ; + + double kFactor = mparams->kFactor(); + + int nbTriangles=_topology->getNbTriangles(); + df.resize(dp.size()); + + for (int i=0; i +void TriangularBendingFEMForceField::convertStiffnessMatrixToGlobalSpace(StiffnessMatrixGlobalSpace &K_gs, TriangleInformation *tinfo) +{ + // Stiffness matrix of current triangle + const StiffnessMatrix &K = tinfo->stiffnessMatrix; + + // Firstly, add all degrees of freedom (we add the unused translation in z) + StiffnessMatrixGlobalSpace K_18x18; + K_18x18.clear(); + unsigned int ig, jg; + + // Copy the stiffness matrix by block 2x2 into 18x18 matrix (the new index of each bloc into global matrix is a combination of 0, 6 and 12 in indices) + for (unsigned int bx=0; bx<3; bx++) + { + // Global row index + ig = 6*bx; + + for (unsigned int by=0; by<3; by++) + { + // Global column index + jg = 6*by; + + // Iterates over the indices of the bloc 2x2 + for (unsigned int i=0; i<2; i++) + { + for (unsigned int j=0; j<2; j++) + { + K_18x18[ig+i][jg+j] = K[2*bx+i][2*by+j]; + } + } + } + } + + + if (d_bending.getValue()) + { + // Stiffness matrix in bending of current triangle + const StiffnessMatrixBending &K_bending = tinfo->stiffnessMatrixBending; + + // Copy the stiffness matrix by block 3x3 into global matrix (the new index of each bloc into global matrix is a combination of 2, 8 and 15 in indices) + for (unsigned int bx=0; bx<3; bx++) + { + // Global row index + ig = 6*bx+2; + + for (unsigned int by=0; by<3; by++) + { + // Global column index + jg = 6*by+2; + + // Iterates over the indices of the bloc 3x3 + for (unsigned int i=0; i<3; i++) + { + for (unsigned int j=0; j<3; j++) + { + K_18x18[ig+i][jg+j] += K_bending[3*bx+i][3*by+j]; + } + } + + } + } + + } + + // Extend rotation matrix and its transpose + Transformation R, Rt; + tinfo->Qframe.toMatrix(R); + Rt.transpose(R); + + StiffnessMatrixGlobalSpace R18x18, Rt18x18; + + for(unsigned int i=0;i<3;++i) + { + for(unsigned int j=0;j<3;++j) + { + R18x18[i][j] = R18x18[i+3][j+3] = R18x18[i+6][j+6] = R18x18[i+9][j+9] = R18x18[i+12][j+12] = R18x18[i+15][j+15] = R[i][j]; + Rt18x18[i][j] = Rt18x18[i+3][j+3] = Rt18x18[i+6][j+6] = Rt18x18[i+9][j+9] = Rt18x18[i+12][j+12] = Rt18x18[i+15][j+15] = Rt[i][j]; + } + } + + // Then we put the stifness matrix into the global frame + K_gs = Rt18x18 * K_18x18 * R18x18; + +} + +template +void TriangularBendingFEMForceField::addKToMatrix(const sofa::core::MechanicalParams* mparams, const sofa::core::behavior::MultiMatrixAccessor* matrix) +{ + StiffnessMatrixGlobalSpace K_gs; + + // Build Matrix Block for this ForceField + unsigned int i, j ,n1, n2, row, column, ROW, COLUMN; + Index node1, node2; + + sofa::core::behavior::MultiMatrixAccessor::MatrixRef r = matrix->getMatrix(this->mstate); + sofa::type::vector& triangleInf = *(triangleInfo.beginEdit()); + + double kFactor = mparams->kFactor(); + + for(sofa::Index t=0 ; t != _topology->getNbTriangles() ; ++t) + { + TriangleInformation *tinfo = &triangleInf[t]; + const Triangle triangle = _topology->getTriangle(t); + + convertStiffnessMatrixToGlobalSpace(K_gs, tinfo); + + // find index of node 1 + for (n1=0; n1<3; n1++) + { + node1 = triangle[n1]; + + for(i=0; i<6; i++) + { + ROW = r.offset+6*node1+i; + row = 6*n1+i; + // find index of node 2 + for (n2=0; n2<3; n2++) + { + node2 = triangle[n2]; + + for (j=0; j<6; j++) + { + COLUMN = r.offset+6*node2+j; + column = 6*n2+j; + r.matrix->add(ROW, COLUMN, - K_gs[row][column] * kFactor); + } + } + } + } + } + + triangleInfo.endEdit(); +} + + +template +void TriangularBendingFEMForceField::addBToMatrix(sofa::linearalgebra::BaseMatrix * /*mat*/, double /*bFact*/, unsigned int &/*offset*/) +{ +} + + +template +void TriangularBendingFEMForceField::draw(const sofa::core::visual::VisualParams* vparams) +{ + if (!vparams->displayFlags().getShowForceFields()) + return; + + const auto stateLifeCycle = vparams->drawTool()->makeStateLifeCycle(); + vparams->drawTool()->disableLighting(); + + if (vparams->displayFlags().getShowWireFrame()) + vparams->drawTool()->setPolygonMode(0, true); + + std::vector colorVector; + std::vector vertices; + + const VecCoord& x = this->mstate->read(sofa::core::vec_id::read_access::position)->getValue(); + + const SeqTriangles triangles = _topology->getTriangles(); + for (const Triangle& tri: triangles) + { + colorVector.push_back(sofa::type::RGBAColor::green()); + vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(x[tri[0]]))); + colorVector.push_back(sofa::type::RGBAColor(0, 0.5, 0.5, 1)); + vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(x[tri[1]]))); + colorVector.push_back(sofa::type::RGBAColor(0, 0, 1, 1)); + vertices.push_back(sofa::type::Vec3(DataTypes::getCPos(x[tri[2]]))); + } + vparams->drawTool()->drawTriangles(vertices, colorVector); + + if (m_mapTopology){ + // Draw lines to visualize the mapping between nodes + std::vector > points; + + // Gets vertices of rest and initial positions respectively + const VecCoord& x0 = (l_restShape.get() != nullptr) + // if having changing rest shape take it + ? l_restShape.get()->f_position.getValue() + : (m_mapTopology + // if rest shape is fixed but we have mapped topology use it + ? l_topologyMapper.get()->f_input_position.getValue() + // otherwise just take rest shape in mechanical state + : this->mstate->read(sofa::core::vec_id::read_access::position)->getValue() + ); + const VecCoord& x = this->mstate->read(sofa::core::vec_id::read_access::position)->getValue(); + + int nbTriangles=_topology->getNbTriangles(); + + for (sofa::Index i=0; i 1e-8) { + points.push_back(x[tinfo.a].getCenter()); + points.push_back(x0[tinfo.a0].getCenter()); + } + if ((x[tinfo.b].getCenter() - x0[tinfo.b0].getCenter()).norm() > 1e-8) { + points.push_back(x[tinfo.b].getCenter()); + points.push_back(x0[tinfo.b0].getCenter()); + } + if ((x[tinfo.c].getCenter() - x0[tinfo.c0].getCenter()).norm() > 1e-8) { + points.push_back(x[tinfo.c].getCenter()); + points.push_back(x0[tinfo.c0].getCenter()); + } + } + + if (points.size() > 0) + vparams->drawTool()->drawLines(points, 1.0f, sofa::type::RGBAColor(1.0, 0.0, 1.0, 1.0)); + } +} + +template +void TriangularBendingFEMForceField::handleEvent(sofa::core::objectmodel::Event *event) +{ + if (dynamic_cast(event)) + { + // Update of the rest shape + // NOTE: the number of triangles should be the same in all topologies + unsigned int nbTriangles = _topology->getNbTriangles(); + for (unsigned int i=0; i. * +******************************************************************************* +* Authors: The SOFA Team and external contributors (see Authors.txt) * +* * +* Contact information: contact@sofa-framework.org * +******************************************************************************/ +#include + +#include +using sofa::core::ObjectFactory; + +extern "C" { + SOFA_SHELL_API void initExternalModule(); + SOFA_SHELL_API const char* getModuleName(); + SOFA_SHELL_API const char* getModuleVersion(); + SOFA_SHELL_API const char* getModuleLicense(); + SOFA_SHELL_API const char* getModuleDescription(); + SOFA_SHELL_API const char* getModuleComponentList(); +} + +void initExternalModule() +{ + static bool first = true; + if (first) + { + first = false; + } +} + +const char* getModuleName() +{ + return sofa_tostring(SOFA_TARGET); +} + +const char* getModuleVersion() +{ + return sofa_tostring(SHELL_VERSION); +} + +const char* getModuleLicense() +{ + return "LGPL"; +} + +const char* getModuleDescription() +{ + return "This plugin implements mechanical models implemented following the Shell method."; +} + +const char* getModuleComponentList() +{ + /// string containing the names of the classes provided by the plugin + static std::string classes = ObjectFactory::getInstance()->listClassesFromTarget(sofa_tostring(SOFA_TARGET)); + return classes.c_str(); +} +