Try the 'lgtm' CodeQL queries

Windows CI METIS workaround

.github/workflows/build-container.yaml

@@ -1,4 +1,4 @@
-name: Build Deploy Container
+name: "Docker"
 
 on:
 

.github/workflows/builds-and-tests.yml

@@ -10,7 +10,7 @@
 # CONTRIBUTING.md for details.
 
 # In this CI section, we build different variants of mfem and run test on them.
-name: builds-and-tests
+name: "Tests"
 
 # Github actions can use the default "GITHUB_TOKEN". By default, this token
 # is set to have permissive access. However, this is not a good practice
@@ -47,17 +47,17 @@ jobs:
   builds-and-tests:
     strategy:
       matrix:
-        os: [ubuntu-20.04, macos-10.15, windows-2022]
+        os: [ubuntu-latest, macos-latest, windows-latest]
         target: [dbg, opt]
         mpi: [seq, par]
         build-system: [make, cmake]
         hypre-target: [int32]
         exclude:
-          - os: ubuntu-20.04
+          - os: ubuntu-latest
             build-system: cmake
-          - os: macos-10.15
+          - os: macos-latest
             build-system: cmake
-          - os: windows-2022
+          - os: windows-latest
             build-system: make
         # 'include' allows us to:
         # - Add a variable to all jobs without creating a new matrix dimension.
@@ -72,15 +72,15 @@ jobs:
             codecov: NO
           - target: opt
             codecov: YES
-          - os: windows-2022
+          - os: windows-latest
             codecov: NO
-          - os: ubuntu-20.04
+          - os: ubuntu-latest
             target: opt
             codecov: NO
             mpi: par
             build-system: cmake
             hypre-target: int32
-          - os: ubuntu-20.04
+          - os: ubuntu-latest
             target: opt
             codecov: NO
             mpi: par
@@ -112,35 +112,35 @@ jobs:
     # TODO: It would be nice to have only one step, e.g. with a dedicated
     # action, but I (@adrienbernede) don't see how at the moment.
     - name: get MPI (Linux)
-      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-20.04'
+      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-latest'
       run: |
         sudo apt-get install mpich libmpich-dev
         export MAKE_CXX_FLAG="MPICXX=mpic++"
 
     - name: get lcov (Linux)
-      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-20.04'
+      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-latest'
       run: |
         sudo apt-get install lcov
 
     - name: Set up Homebrew
-      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-10.15'
+      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
       uses: Homebrew/actions/setup-homebrew@c4aafe8c4620bf08883dd4679c374f11e73329d3
 
     - name: get MPI (MacOS)
-      if: matrix.mpi == 'par' && matrix.os == 'macos-10.15'
+      if: matrix.mpi == 'par' && matrix.os == 'macos-latest'
       run: |
         export HOMEBREW_NO_INSTALL_CLEANUP=1
         brew install openmpi
         export MAKE_CXX_FLAG="MPICXX=mpic++"
 
     - name: get MPI (MacOS)
-      if: matrix.codecov == 'YES' && matrix.os == 'macos-10.15'
+      if: matrix.codecov == 'YES' && matrix.os == 'macos-latest'
       run: |
         export HOMEBREW_NO_INSTALL_CLEANUP=1
         brew install lcov
 
     - name: get MPI (Windows)
-      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
+      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
       uses: mpi4py/setup-mpi@v1.0.3
 
     # Get Hypre through cache, or build it.
@@ -154,7 +154,7 @@ jobs:
         key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.2
 
     - name: get hypre
-      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-2022'
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-latest'
       uses: mfem/github-actions/build-hypre@v2.2
       with:
         archive: ${{ env.HYPRE_ARCHIVE }}
@@ -163,7 +163,7 @@ jobs:
         build-system: make
 
     - name: get hypre (Windows)
-      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-2022'
+      if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-latest'
       uses: mfem/github-actions/build-hypre@v2.2
       with:
         archive: ${{ env.HYPRE_ARCHIVE }}
@@ -175,14 +175,14 @@ jobs:
     # Install will only run on cache miss.
     - name: cache metis
       id: metis-cache
-      if: matrix.mpi == 'par' && matrix.os != 'windows-2022'
+      if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
       uses: actions/cache@v2
       with:
         path: ${{ env.METIS_TOP_DIR }}
         key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
 
     - name: install metis
-      if: matrix.mpi == 'par' && matrix.os != 'windows-2022' && steps.metis-cache.outputs.cache-hit != 'true'
+      if: matrix.mpi == 'par' && matrix.os != 'windows-latest' && steps.metis-cache.outputs.cache-hit != 'true'
       uses: mfem/github-actions/build-metis@v2.2
       with:
         archive: ${{ env.METIS_ARCHIVE }}
@@ -196,20 +196,16 @@ jobs:
         key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
 
     - name: prepare binary cache location
-      if: matrix.os == 'windows-2022' && steps.vcpkg-cache.outputs.cache-hit != 'true'
+      if: matrix.os == 'windows-latest' && steps.vcpkg-cache.outputs.cache-hit != 'true'
       run: |
         mkdir -p vcpkg_cache
 
     - name: install metis (Windows)
-      if: matrix.mpi == 'par' && matrix.os == 'windows-2022'
+      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
       env:
         VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
       run: |
-        $PortFile = 'C:\vcpkg\ports\metis\portfile.cmake'
-        $OriginalURL = 'http://glaros.dtc.umn.edu/gkhome/fetch/sw/metis/metis-${METIS_VERSION}.tar.gz'
-        $NewURL = 'https://github.com/mfem/tpls/raw/gh-pages/metis-5.1.0.tar.gz'
-        (Get-Content $PortFile).replace($OriginalURL, $NewURL) | Set-Content $PortFile
-        vcpkg install metis --triplet=x64-windows-static
+        vcpkg install metis-mfem --triplet=x64-windows-static --overlay-ports=${{ env.MFEM_TOP_DIR }}/config/vcpkg/ports
 
     # MFEM build and test
     - name: build
@@ -252,7 +248,7 @@ jobs:
       shell: bash
 
     - name: cmake unit tests (Ubuntu 20.04)
-      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-20.04'
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-latest'
       run: |
         CTEST_CONFIG="Release"
         [[ ${{ matrix.target }} == 'dbg' ]] && CTEST_CONFIG="Debug"
@@ -260,7 +256,7 @@ jobs:
       shell: bash
 
     - name: cmake tests
-      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-20.04'
+      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os != 'ubuntu-latest'
       run: |
         CTEST_CONFIG="Release"
         cd ${{ env.MFEM_TOP_DIR }}/build && ctest --output-on-failure -C ${CTEST_CONFIG}

.github/workflows/codeql-analysis.yml

@@ -0,0 +1,71 @@
+# For most projects, this workflow file will not need changing; you simply need
+# to commit it to your repository.
+#
+# You may wish to alter this file to override the set of languages analyzed,
+# or to provide custom queries or build logic.
+#
+# ******** NOTE ********
+# We have attempted to detect the languages in your repository. Please check
+# the `language` matrix defined below to confirm you have the correct set of
+# supported CodeQL languages.
+#
+name: "Static Analysis"
+
+on:
+  push:
+    branches: [ "master", "next"]
+  pull_request:
+    # The branches below must be a subset of the branches above
+    branches: [ "master" ]
+
+jobs:
+  analyze:
+    name: Analyze
+    runs-on: ubuntu-latest
+    permissions:
+      actions: read
+      contents: read
+      security-events: write
+
+    strategy:
+      fail-fast: false
+      matrix:
+        language: [ 'cpp' ]
+        # CodeQL supports [ 'cpp', 'csharp', 'go', 'java', 'javascript', 'python', 'ruby' ]
+        # Learn more about CodeQL language support at https://aka.ms/codeql-docs/language-support
+
+    steps:
+    - name: Checkout repository
+      uses: actions/checkout@v3
+
+    # Initializes the CodeQL tools for scanning.
+    - name: Initialize CodeQL
+      uses: github/codeql-action/init@v2
+      with:
+        languages: ${{ matrix.language }}
+        # If you wish to specify custom queries, you can do so here or in a config file.
+        # By default, queries listed here will override any specified in a config file.
+        # Prefix the list here with "+" to use these queries and those in the config file.
+
+        # Details on CodeQL's query packs refer to : https://docs.github.com/en/code-security/code-scanning/automatically-scanning-your-code-for-vulnerabilities-and-errors/configuring-code-scanning#using-queries-in-ql-packs
+        # queries: security-extended,security-and-quality
+        queries: lgtm
+
+
+    # Autobuild attempts to build any compiled languages  (C/C++, C#, or Java).
+    # If this step fails, then you should remove it and run the build manually (see below)
+    - name: Autobuild
+      uses: github/codeql-action/autobuild@v2
+
+    # ℹ️ Command-line programs to run using the OS shell.
+    # 📚 See https://docs.github.com/en/actions/using-workflows/workflow-syntax-for-github-actions#jobsjob_idstepsrun
+
+    #   If the Autobuild fails above, remove it and uncomment the following three lines.
+    #   modify them (or add more) to build your code if your project, please refer to the EXAMPLE below for guidance.
+
+    # - run: |
+    #   echo "Run, Build Application using script"
+    #   ./location_of_script_within_repo/buildscript.sh
+
+    - name: Perform CodeQL Analysis
+      uses: github/codeql-action/analyze@v2

.github/workflows/mfem-analysis.yml

@@ -9,7 +9,7 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-name: build-analysis
+name: "Build Analysis"
 
 permissions:
   actions: write
@@ -31,7 +31,7 @@ env:
 
 jobs:
   gitignore:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
 
     steps:
     - name: Cancel Previous Runs

.github/workflows/repo-check.yml

@@ -9,7 +9,7 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-name: repo-check
+name: "Checks"
 
 permissions:
   actions: write
@@ -28,7 +28,7 @@ on:
 
 jobs:
   file-headers-check:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -78,7 +78,7 @@ jobs:
         exit 1
 
   code-style:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -88,14 +88,14 @@ jobs:
 
     - name: get astyle
       run: |
-        sudo apt-get install astyle=3.1-1ubuntu2
+        sudo apt-get install astyle
 
     - name: style check
       run: |
         ./config/githooks/pre-push --style
 
   documentation:
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     if: |
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
@@ -106,6 +106,8 @@ jobs:
     - name: get doxygen and graphviz
       run: |
         sudo apt-get install doxygen graphviz
+        cd doc
+        doxygen -u CodeDocumentation.conf.in 2>/dev/null
 
     - name: build documentation
       run: |
@@ -118,7 +120,7 @@ jobs:
       github.ref != 'refs/heads/master' &&
       (github.event_name == 'push' ||
        github.event.pull_request.head.repo.full_name != github.repository)
-    runs-on: ubuntu-18.04
+    runs-on: ubuntu-latest
     steps:
     - name: checkout mfem
       uses: actions/checkout@v2

config/vcpkg/ports/metis-mfem/disable-programs.patch

@@ -0,0 +1,8 @@
+--- a/CMakeLists.txt	Wed Dec 21 18:24:22 2016
++++ b/CMakeLists.txt	Wed Dec 21 18:24:26 2016
+@@ -20,4 +20,4 @@
+ # Recursively look for CMakeLists.txt in subdirs.
+ add_subdirectory("include")
+ add_subdirectory("libmetis")
+-add_subdirectory("programs")
++# add_subdirectory("programs")

config/vcpkg/ports/metis-mfem/enable-install.patch

@@ -0,0 +1,15 @@
+--- a/CMakeLists.txt	Sat Mar 30 17:24:45 2013
++++ b/CMakeLists.txt	Wed Dec 21 18:23:43 2016
+@@ -4,11 +4,7 @@
+ set(GKLIB_PATH "GKlib" CACHE PATH "path to GKlib")
+ set(SHARED FALSE CACHE BOOL "build a shared library")
+ 
+-if(MSVC)
+-  set(METIS_INSTALL FALSE)
+-else()
+-  set(METIS_INSTALL TRUE)
+-endif()
++set(METIS_INSTALL TRUE)
+ 
+ # Configure libmetis library.
+ if(SHARED)

config/vcpkg/ports/metis-mfem/fix-INT_MIN_define.patch

@@ -0,0 +1,34 @@
+diff --git a/include/metis.h b/include/metis.h
+index dc5406a..7732437 100644
+--- a/include/metis.h
++++ b/include/metis.h
+@@ -72,10 +72,14 @@ typedef __int64 int64_t;
+ #define PRId64       "I64d"
+ #define SCNd32       "ld"
+ #define SCNd64       "I64d"
++#ifdef _WIN32
++#include <stdint.h>
++#else
+ #define INT32_MIN    ((int32_t)_I32_MIN)
+ #define INT32_MAX    _I32_MAX
+ #define INT64_MIN    ((int64_t)_I64_MIN)
+ #define INT64_MAX    _I64_MAX
++#endif
+ #else
+ #include <inttypes.h>
+ #endif
+diff --git a/GKlib/gk_arch.h b/GKlib/gk_arch.h
+index 78b1431..7258763 100644
+--- a/GKlib/gk_arch.h
++++ b/GKlib/gk_arch.h
+@@ -32,8 +32,8 @@
+ 
+ 
+ #ifdef __MSC__ 
+-  #include "ms_stdint.h"
+-  #include "ms_inttypes.h"
++  #include <stdint.h>
++  #include <inttypes.h>
+   #include "ms_stat.h"
+ #else
+ #ifndef SUNOS

config/vcpkg/ports/metis-mfem/fix-gklib-vs14-math.patch

@@ -0,0 +1,11 @@
+--- a/GKlib/gk_arch.h	Wed Dec 21 18:34:18 2016
++++ b/GKlib/gk_arch.h	Wed Dec 21 18:30:49 2016
+@@ -58,7 +58,7 @@
+ #define PTRDIFF_MAX  INT64_MAX
+ #endif
+ 
+-#ifdef __MSC__
++#if defined(__MSC__) && (_MSC_VER < 1900)
+ /* MSC does not have rint() function */
+ #define rint(x) ((int)((x)+0.5))  
+ 

config/vcpkg/ports/metis-mfem/fix-linux-build-error.patch

@@ -0,0 +1,14 @@
+diff --git a/CMakeLists.txt b/CMakeLists.txt
+index e94f050..b9613a7 100644
+--- a/CMakeLists.txt
++++ b/CMakeLists.txt
+@@ -1,7 +1,8 @@
+ cmake_minimum_required(VERSION 2.8)
+ project(METIS)
+ 
+-set(GKLIB_PATH "GKlib" CACHE PATH "path to GKlib")
++set(GKLIB_PATH "${CMAKE_SOURCE_DIR}/GKlib" CACHE PATH "path to GKlib")
++
+ set(SHARED FALSE CACHE BOOL "build a shared library")
+ 
+ set(METIS_INSTALL TRUE)

config/vcpkg/ports/metis-mfem/fix-metis-vs14-math.patch

@@ -0,0 +1,11 @@
+--- a/libmetis/metislib.h	Sat Mar 30 17:24:45 2013
++++ b/libmetis/metislib.h	Wed Dec 21 18:30:59 2016
+@@ -31,7 +31,7 @@
+ #include <proto.h>
+ 
+ 
+-#if defined(COMPILER_MSC)
++#if defined(COMPILER_MSC) && (_MSC_VER < 1900)
+ #if defined(rint)
+   #undef rint
+ #endif

config/vcpkg/ports/metis-mfem/fix-runtime-install-destination.patch

@@ -0,0 +1,10 @@
+--- a/libmetis/CMakeLists.txt	Sat Mar 30 17:24:45 2013
++++ b/libmetis/CMakeLists.txt	Wed Dec 21 17:41:37 2016
+@@ -11,6 +11,6 @@
+ if(METIS_INSTALL)
+   install(TARGETS metis
+     LIBRARY DESTINATION lib
+-    RUNTIME DESTINATION lib
++    RUNTIME DESTINATION bin
+     ARCHIVE DESTINATION lib)
+ endif()

config/vcpkg/ports/metis-mfem/install-metisConfig.patch

@@ -0,0 +1,44 @@
+diff --git a/CMakeLists.txt b/CMakeLists.txt
+index b9613a7..e43ffee 100644
+--- a/CMakeLists.txt
++++ b/CMakeLists.txt
+@@ -22,3 +22,23 @@ include_directories(include)
+ add_subdirectory("include")
+ add_subdirectory("libmetis")
+ # add_subdirectory("programs")
++
++if(METIS_INSTALL)
++  set(PRJ_NAME metis)
++  set(PRJ_VER 5.1.0)
++  install(EXPORT metisTargets 
++    FILE ${PRJ_NAME}Targets.cmake
++    DESTINATION lib/cmake/${PRJ_NAME})
++  include(CMakePackageConfigHelpers)
++  write_basic_package_version_file(
++    ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}ConfigVersion.cmake
++    VERSION ${PRJ_VER}
++    COMPATIBILITY SameMajorVersion)
++  file(WRITE ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}Config.cmake
++    "include(\${CMAKE_CURRENT_LIST_DIR}/${PRJ_NAME}Targets.cmake)")
++  install(FILES 
++    ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}ConfigVersion.cmake
++    ${CMAKE_CURRENT_BINARY_DIR}/${PRJ_NAME}Config.cmake
++    DESTINATION lib/cmake/${PRJ_NAME})
++endif()
++
+diff --git a/libmetis/CMakeLists.txt b/libmetis/CMakeLists.txt
+index 7a5fc74..5a68cf0 100644
+--- a/libmetis/CMakeLists.txt
++++ b/libmetis/CMakeLists.txt
+@@ -9,8 +9,9 @@ if(UNIX)
+ endif()
+ 
+ if(METIS_INSTALL)
+-  install(TARGETS metis
++  install(TARGETS metis EXPORT metisTargets
+     LIBRARY DESTINATION lib
+     RUNTIME DESTINATION bin
+-    ARCHIVE DESTINATION lib)
++    ARCHIVE DESTINATION lib
++    INCLUDES DESTINATION include)
+ endif()

config/vcpkg/ports/metis-mfem/portfile.cmake

@@ -0,0 +1,41 @@
+vcpkg_check_linkage(ONLY_STATIC_LIBRARY)
+set(OPTIONS -DSHARED=OFF)
+
+set(METIS_VERSION 5.1.0)
+
+vcpkg_download_distfile(ARCHIVE
+    URLS "https://github.com/mfem/tpls/raw/gh-pages/metis-${METIS_VERSION}.tar.gz"
+    FILENAME "metis-${METIS_VERSION}.tar.gz"
+    SHA512 deea47749d13bd06fbeaf98a53c6c0b61603ddc17a43dae81d72c8015576f6495fd83c11b0ef68d024879ed5415c14ebdbd87ce49c181bdac680573bea8bdb25
+)
+
+vcpkg_extract_source_archive_ex(
+    OUT_SOURCE_PATH SOURCE_PATH
+    ARCHIVE ${ARCHIVE}
+    REF ${METIS_VERSION}
+    PATCHES
+        enable-install.patch
+        disable-programs.patch
+        fix-runtime-install-destination.patch
+        fix-metis-vs14-math.patch
+        fix-gklib-vs14-math.patch
+        fix-linux-build-error.patch
+        install-metisConfig.patch
+        fix-INT_MIN_define.patch
+)
+
+vcpkg_configure_cmake(
+    SOURCE_PATH ${SOURCE_PATH}
+    PREFER_NINJA
+    OPTIONS ${OPTIONS}
+)
+
+vcpkg_install_cmake()
+vcpkg_copy_pdbs()
+vcpkg_fixup_cmake_targets(CONFIG_PATH lib/cmake/metis)
+
+file(REMOVE_RECURSE ${CURRENT_PACKAGES_DIR}/debug/include)
+
+# Handle copyright
+file(COPY ${SOURCE_PATH}/LICENSE.txt DESTINATION ${CURRENT_PACKAGES_DIR}/share/metis)
+file(INSTALL ${SOURCE_PATH}/LICENSE.txt DESTINATION ${CURRENT_PACKAGES_DIR}/share/${PORT} RENAME copyright)

config/vcpkg/ports/metis-mfem/vcpkg.json

@@ -0,0 +1,7 @@
+{
+  "name": "metis-mfem",
+  "version-string": "5.1.0",
+  "port-version": 0,
+  "description": "Serial Graph Partitioning and Fill-reducing Matrix Ordering",
+  "homepage": "https://glaros.dtc.umn.edu/gkhome/metis/metis/overview"
+}

doc/CodeDocumentation.conf.in

@@ -2349,7 +2349,7 @@ PLANTUML_INCLUDE_PATH  =
 # Minimum value: 0, maximum value: 10000, default value: 50.
 # This tag requires that the tag HAVE_DOT is set to YES.
 
-DOT_GRAPH_MAX_NODES    = 50
+DOT_GRAPH_MAX_NODES    = 100
 
 # The MAX_DOT_GRAPH_DEPTH tag can be used to set the maximum depth of the graphs
 # generated by dot. A depth value of 3 means that only nodes reachable from the

examples/ex1.cpp

@@ -30,7 +30,7 @@
 //
 // Device sample runs:
 //               ex1 -pa -d cuda
-//               ex1 -fa -d cuda
+//             * ex1 -fa -d cuda
 //               ex1 -pa -d raja-cuda
 //             * ex1 -pa -d raja-hip
 //               ex1 -pa -d occa-cuda

examples/ex1p.cpp

@@ -30,7 +30,7 @@
 //
 // Device sample runs:
 //               mpirun -np 4 ex1p -pa -d cuda
-//               mpirun -np 4 ex1p -fa -d cuda
+//             * mpirun -np 4 ex1p -fa -d cuda
 //               mpirun -np 4 ex1p -pa -d occa-cuda
 //               mpirun -np 4 ex1p -pa -d raja-omp
 //               mpirun -np 4 ex1p -pa -d ceed-cpu

fem/lor/lor.cpp

@@ -467,7 +467,8 @@ void LORDiscretization::FormLORSpace()
    mesh = new Mesh(Mesh::MakeRefined(mesh_ho, refinements, ref_type));
 
    fec = fes_ho.FEColl()->Clone(GetLOROrder());
-   fes = new FiniteElementSpace(mesh, fec);
+   const int vdim = fes_ho.GetVDim();
+   fes = new FiniteElementSpace(mesh, fec, vdim);
    SetupProlongationAndRestriction();
 }
 
@@ -511,7 +512,8 @@ void ParLORDiscretization::FormLORSpace()
    mesh = pmesh;
 
    fec = pfes_ho.FEColl()->Clone(GetLOROrder());
-   ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec);
+   const int vdim = fes_ho.GetVDim();
+   ParFiniteElementSpace *pfes = new ParFiniteElementSpace(pmesh, fec, vdim);
    fes = pfes;
    SetupProlongationAndRestriction();
 }

fem/tbilininteg.hpp

@@ -361,8 +361,8 @@ struct TDiffusionKernel<2,2,complex_t>
        @param F Jt [M x Dim x SDim x NE] - Jacobian transposed, data member in F
        @param Q CoefficientEval<>::Type
        @param q CoefficientEval<>::Type::result_t
-       @param A [M x Dim*(Dim+1)/2] partially assembled Dim x Dim symm. matrices
-       @param A [M x Dim x Dim]       partially assembled Dim x Dim matrices
+       @param A either [M x Dim*(Dim+1)/2] partially assembled Dim x Dim symm.
+       matrices, or [M x Dim x Dim] partially assembled Dim x Dim matrices.
    */
    template <typename T_result_t, typename Q_t, typename q_t, typename asm_type>
    static inline MFEM_ALWAYS_INLINE

linalg/ode.cpp

@@ -188,7 +188,6 @@ void ExplicitRKSolver::Step(Vector &x, double &t, double &dt)
 
 ExplicitRKSolver::~ExplicitRKSolver()
 {
-   delete [] k;
 }
 
 const double RK6Solver::a[] =

linalg/ode.hpp

@@ -331,7 +331,6 @@ public:
    ~AdamsMoultonSolver()
    {
       if (RKsolver) { delete RKsolver; }
-      delete [] k;
    };
 };
 

miniapps/CMakeLists.txt

@@ -31,5 +31,4 @@ add_subdirectory(shifted)
 add_subdirectory(mtop)
 add_subdirectory(autodiff)
 add_subdirectory(parelag)
-add_subdirectory(elplast)
 add_subdirectory(hooke)

miniapps/elplast/CMakeLists.txt

@@ -1,49 +0,0 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-if (MFEM_USE_MPI)
-  list(APPEND DIST_COMMON_SOURCES
-	  coefficients.cpp
-	  elplast.cpp
-      )
-  list(APPEND DIST_COMMON_HEADERS
-	  coefficients.hpp
-	  elplast.hpp
-      )
-
-  convert_filenames_to_full_paths(DIST_COMMON_SOURCES)
-  convert_filenames_to_full_paths(DIST_COMMON_HEADERS)
-
-  set(DIST_COMMON_FILES
-      EXTRA_SOURCES ${DIST_COMMON_SOURCES}
-      EXTRA_HEADERS ${DIST_COMMON_HEADERS})
-
-  add_mfem_miniapp(adapt
-    MAIN adaptive_el.cpp
-    ${DIST_COMMON_FILES}
-    LIBRARIES mfem mfem-common)
-
-  add_mfem_miniapp(test_elplinteg
-    MAIN test_elinteg.cpp
-    ${DIST_COMMON_FILES}
-    LIBRARIES mfem mfem-common)
-
-  add_mfem_miniapp(elpl_elem
-    MAIN elpl_test.cpp
-    ${DIST_COMMON_FILES}
-    LIBRARIES mfem mfem-common)
-
-  add_mfem_miniapp(holep
-    MAIN holep.cpp
-    ${DIST_COMMON_FILES}
-    LIBRARIES mfem mfem-common)
-
-endif ()

miniapps/elplast/adaptive_el.cpp

@@ -1,401 +0,0 @@
-//                       MFEM Example 21 - Parallel Version
-//
-// Compile with: make ex21p
-//
-// Sample runs:  mpirun -np 4 ex21p
-//               mpirun -np 4 ex21p -o 3
-//               mpirun -np 4 ex21p -m ../data/beam-quad.mesh
-//               mpirun -np 4 ex21p -m ../data/beam-quad.mesh -o 3
-//               mpirun -np 4 ex21p -m ../data/beam-tet.mesh
-//               mpirun -np 4 ex21p -m ../data/beam-tet.mesh -o 2
-//               mpirun -np 4 ex21p -m ../data/beam-hex.mesh
-//               mpirun -np 4 ex21p -m ../data/beam-hex.mesh -o 2
-//
-// Description:  This is a version of Example 2p with a simple adaptive mesh
-//               refinement loop. The problem being solved is again the linear
-//               elasticity describing a multi-material cantilever beam.
-//               The problem is solved on a sequence of meshes which
-//               are locally refined in a conforming (triangles, tetrahedrons)
-//               or non-conforming (quadrilaterals, hexahedra) manner according
-//               to a simple ZZ error estimator.
-//
-//               The example demonstrates MFEM's capability to work with both
-//               conforming and nonconforming refinements, in 2D and 3D, on
-//               linear and curved meshes. Interpolation of functions from
-//               coarse to fine meshes, as well as persistent GLVis
-//               visualization are also illustrated.
-//
-//               We recommend viewing Examples 2p and 6p before viewing this
-//               example.
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-#include "coefficients.hpp"
-
-using namespace std;
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   // 0. Initialize MPI and HYPRE.
-   Mpi::Init(argc, argv);
-   int num_procs = Mpi::WorldSize();
-   int myid = Mpi::WorldRank();
-   Hypre::Init();
-
-   // 1. Parse command-line options.
-   const char *mesh_file = "../data/beam-tri.mesh";
-   int serial_ref_levels = 0;
-   int order = 2;
-   //bool static_cond = false;
-   bool static_cond = true;
-   bool visualization = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&serial_ref_levels, "-rs", "--refine-serial",
-                  "Number of uniform serial refinements (before parallel"
-                  " partitioning)");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree).");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   // 2. Read the mesh from the given mesh file. We can handle triangular,
-   //    quadrilateral, tetrahedral, and hexahedral meshes with the same code.
-   Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-   MFEM_VERIFY(mesh.SpaceDimension() == dim, "invalid mesh");
-
-   // 3. Refine the mesh before parallel partitioning. Since a NURBS mesh can
-   //    currently only be refined uniformly, we need to convert it to a
-   //    piecewise-polynomial curved mesh. First we refine the NURBS mesh a bit
-   //    more and then project the curvature to quadratic Nodes.
-   if (mesh.NURBSext && serial_ref_levels == 0)
-   {
-      serial_ref_levels = 2;
-   }
-   for (int i = 0; i < serial_ref_levels; i++)
-   {
-      mesh.UniformRefinement();
-   }
-   if (mesh.NURBSext)
-   {
-      mesh.SetCurvature(2);
-   }
-   mesh.EnsureNCMesh();
-
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   mesh.Clear();
-
-   // 4. Define a finite element space on the mesh. The polynomial order is
-   //    one (linear) by default, but this can be changed on the command line.
-   H1_FECollection fec(order, dim);
-   ParFiniteElementSpace fespace(&pmesh, &fec, dim);
-
-   // 5. As in Example 2, we set up the linear form b(.) which corresponds to
-   //    the right-hand side of the FEM linear system. In this case, b_i equals
-   //    the boundary integral of f*phi_i where f represents a "pull down"
-   //    force on the Neumann part of the boundary and phi_i are the basis
-   //    functions in the finite element fespace. The force is defined by the
-   //    VectorArrayCoefficient object f, which is a vector of Coefficient
-   //    objects. The fact that f is non-zero on boundary attribute 2 is
-   //    indicated by the use of piece-wise constants coefficient for its last
-   //    component. We don't assemble the discrete problem yet, this will be
-   //    done in the main loop.
-   VectorArrayCoefficient f(dim);
-   f.Set(0, new ConstantCoefficient(0.0));
-   f.Set(1, new ConstantCoefficient(1.0));
-   if(dim==3){
-       f.Set(2, new ConstantCoefficient(0.0));
-   }
-
-   ParLinearForm b(&fespace);
-
-   b.AddDomainIntegrator(new VectorDomainLFIntegrator(f));
-
-
-   // 6. Set up the bilinear form a(.,.) on the finite element space
-   //    corresponding to the linear elasticity integrator with piece-wise
-   //    constants coefficient lambda and mu.
-   /*
-   Vector lambda(pmesh.attributes.Max());
-   lambda = 1.0;
-   lambda(0) = lambda(1)*50;
-   PWConstCoefficient lambda_func(lambda);
-   Vector mu(pmesh.attributes.Max());
-   mu = 1.0;
-   mu(0) = mu(1)*50;
-   PWConstCoefficient mu_func(mu);
-   */
-
-   //ConstantCoefficient lambda_func(1.0);
-   ConstantCoefficient lambda_func(0.576923E+7);
-   //ConstantCoefficient mu_func(1.0);
-   ConstantCoefficient mu_func(0.384615E+7);
-
-   ParBilinearForm a(&fespace);
-   BilinearFormIntegrator *integ =
-      new ElasticityIntegrator(lambda_func,mu_func);
-   a.AddDomainIntegrator(integ);
-   if (static_cond) { a.EnableStaticCondensation(); }
-
-   // 7. The solution vector x and the associated finite element grid function
-   //    will be maintained over the AMR iterations. We initialize it to zero.
-   Vector zero_vec(dim);
-   zero_vec = 0.0;
-   VectorConstantCoefficient zero_vec_coeff(zero_vec);
-   ParGridFunction x(&fespace);
-   x = 0.0;
-
-   // 8. Determine the list of true (i.e. conforming) essential boundary dofs.
-   //    In this example, the boundary conditions are defined by marking only
-   //    boundary attribute 1 from the mesh as essential and converting it to a
-   //    list of true dofs.  The conversion to true dofs will be done in the
-   //    main loop.
-   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
-   ess_bdr = 0;
-   ess_bdr[0] = 1;
-
-   // 9. GLVis visualization.
-   char vishost[] = "localhost";
-   int  visport   = 19916;
-   socketstream sol_sock;
-
-   // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
-   //     that uses the ComputeElementFlux method of the ElasticityIntegrator to
-   //     recover a smoothed flux (stress) that is subtracted from the element
-   //     flux to get an error indicator. We need to supply the space for the
-   //     smoothed flux: an (H1)^tdim (i.e., vector-valued) space is used here.
-   //     Here, tdim represents the number of components for a symmetric (dim x
-   //     dim) tensor.
-   const int tdim = dim*(dim+1)/2;
-   L2_FECollection flux_fec(order, dim);
-   ParFiniteElementSpace flux_fespace(&pmesh, &flux_fec, tdim);
-   ParFiniteElementSpace smooth_flux_fespace(&pmesh, &fec, tdim);
-   L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace,
-                                       smooth_flux_fespace);
-
-
-   //strain fields
-   L2_FECollection strain_fec(order-1,dim);
-   ParFiniteElementSpace strain_fespace(&pmesh,&strain_fec, dim*(dim+1)/2);
-   ParGridFunction estrains; estrains.SetSpace(&strain_fespace);
-   ParGridFunction sstrains; sstrains.SetSpace(&strain_fespace);
-
-
-   // ParaView output.
-   ParaViewDataCollection dacol("ParaView", &pmesh);
-   dacol.SetLevelsOfDetail(order);
-   dacol.SetDataFormat(VTKFormat::ASCII);
-   dacol.RegisterField("disp", &x);
-   dacol.RegisterField("estr",&estrains);
-   dacol.RegisterField("sstr",&sstrains);
-   dacol.SetTime(1.0);
-   dacol.SetCycle(1);
-   dacol.Save();
-
-   // 11. A refiner selects and refines elements based on a refinement strategy.
-   //     The strategy here is to refine elements with errors larger than a
-   //     fraction of the maximum element error. Other strategies are possible.
-   //     The refiner will call the given error estimator.
-   ThresholdRefiner refiner(estimator);
-   refiner.SetTotalErrorFraction(0.7);
-
-   // 12. The main AMR loop. In each iteration we solve the problem on the
-   //     current mesh, visualize the solution, and refine the mesh.
-   const int max_dofs = 100000;
-   const int max_amr_itr = 20;
-   for (int it = 0; it <= max_amr_itr; it++)
-   {
-      HYPRE_BigInt global_dofs = fespace.GlobalTrueVSize();
-      if (myid == 0)
-      {
-         cout << "\nAMR iteration " << it << endl;
-         cout << "Number of unknowns: " << global_dofs << endl;
-      }
-
-      // 13. Assemble the stiffness matrix and the right-hand side.
-      a.Assemble();
-      b.Assemble();
-
-      // 14. Set Dirichlet boundary values in the GridFunction x.
-      //     Determine the list of Dirichlet true DOFs in the linear system.
-      Array<int> ess_tdof_list;
-      x.ProjectBdrCoefficient(zero_vec_coeff, ess_bdr);
-      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-
-      // 15. Create the linear system: eliminate boundary conditions, constrain
-      //     hanging nodes and possibly apply other transformations. The system
-      //     will be solved for true (unconstrained) DOFs only.
-
-      HypreParMatrix A;
-      Vector B, X;
-      const int copy_interior = 1;
-      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
-
-      // 16. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
-      //     preconditioner from hypre.
-      HypreBoomerAMG amg;
-      amg.SetPrintLevel(0);
-      // amg.SetSystemsOptions(dim); // optional
-      CGSolver pcg(A.GetComm());
-      pcg.SetPreconditioner(amg);
-      pcg.SetOperator(A);
-      pcg.SetRelTol(1e-6);
-      pcg.SetMaxIter(500);
-      pcg.SetPrintLevel(3); // print the first and the last iterations only
-      pcg.Mult(B, X);
-
-      // 17. After solving the linear system, reconstruct the solution as a
-      //     finite element GridFunction. Constrained nodes are interpolated
-      //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
-      a.RecoverFEMSolution(X, b, x);
-
-      //compute the starins
-      EngStrainCoefficient estr(x);
-      StrainCoefficient sstr(x);
-
-      estrains.ProjectCoefficient(estr);
-      sstrains.ProjectCoefficient(sstr);
-
-      //ParaView output
-      dacol.SetTime(double(it+1));
-      dacol.SetCycle(it+1);
-      dacol.Save();
-
-
-      // 18. Send solution by socket to the GLVis server.
-      if (visualization && it == 0)
-      {
-         sol_sock.open(vishost, visport);
-         sol_sock.precision(8);
-      }
-      if (visualization && sol_sock.good())
-      {
-         GridFunction nodes(&fespace), *nodes_p = &nodes;
-         pmesh.GetNodes(nodes);
-         nodes += x;
-         int own_nodes = 0;
-         pmesh.SwapNodes(nodes_p, own_nodes);
-         x.Neg(); // visualize the backward displacement
-         sol_sock << "parallel " << num_procs << ' ' << myid << '\n';
-         sol_sock << "solution\n" << pmesh << x << flush;
-         x.Neg();
-         pmesh.SwapNodes(nodes_p, own_nodes);
-         if (it == 0)
-         {
-            sol_sock << "keys '" << ((dim == 2) ? "Rjl" : "") << "m'" << endl;
-         }
-         sol_sock << "window_title 'AMR iteration: " << it << "'\n"
-                  << "pause" << endl;
-         if (myid == 0)
-         {
-            cout << "Visualization paused. "
-                 "Press <space> in the GLVis window to continue." << endl;
-         }
-      }
-
-      if (global_dofs > max_dofs)
-      {
-         if (myid == 0)
-         {
-            cout << "Reached the maximum number of dofs. Stop." << endl;
-         }
-         break;
-      }
-
-      // 19. Call the refiner to modify the mesh. The refiner calls the error
-      //     estimator to obtain element errors, then it selects elements to be
-      //     refined and finally it modifies the mesh. The Stop() method can be
-      //     used to determine if a stopping criterion was met.
-      refiner.Apply(pmesh);
-      if (refiner.Stop())
-      {
-         if (myid == 0)
-         {
-            cout << "Stopping criterion satisfied. Stop." << endl;
-         }
-         break;
-      }
-
-      // 20. Update the space to reflect the new state of the mesh. Also,
-      //     interpolate the solution x so that it lies in the new space but
-      //     represents the same function. This saves solver iterations later
-      //     since we'll have a good initial guess of x in the next step.
-      //     Internally, FiniteElementSpace::Update() calculates an
-      //     interpolation matrix which is then used by GridFunction::Update().
-      fespace.Update();
-      x.Update();
-
-      strain_fespace.Update();
-      estrains.Update();
-      sstrains.Update();
-
-      // 21. Load balance the mesh, and update the space and solution. Currently
-      //     available only for nonconforming meshes.
-      if (pmesh.Nonconforming())
-      {
-         pmesh.Rebalance();
-
-         // Update the space and the GridFunction. This time the update matrix
-         // redistributes the GridFunction among the processors.
-         fespace.Update();
-         x.Update();
-      }
-
-      // 21. Inform also the bilinear and linear forms that the space has
-      //     changed.
-      a.Update();
-      b.Update();
-   }
-
-   dacol.SetTime(double(max_amr_itr));
-   dacol.SetCycle(max_amr_itr);
-   dacol.Save();
-
-   {
-      ostringstream mref_name, mesh_name, sol_name;
-      mref_name << "ex21p_reference_mesh." << setfill('0') << setw(6) << myid;
-      mesh_name << "ex21p_deformed_mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "ex21p_displacement." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ref_out(mref_name.str().c_str());
-      mesh_ref_out.precision(16);
-      pmesh.Print(mesh_ref_out);
-
-      ofstream mesh_out(mesh_name.str().c_str());
-      mesh_out.precision(16);
-      GridFunction nodes(&fespace), *nodes_p = &nodes;
-      pmesh.GetNodes(nodes);
-      nodes += x;
-      int own_nodes = 0;
-      pmesh.SwapNodes(nodes_p, own_nodes);
-      pmesh.Print(mesh_out);
-      pmesh.SwapNodes(nodes_p, own_nodes);
-
-      ofstream x_out(sol_name.str().c_str());
-      x_out.precision(16);
-      x.Save(x_out);
-   }
-
-   return 0;
-}

miniapps/elplast/coefficients.hpp

@@ -1,733 +0,0 @@
-#ifndef COEFFICIENTS_HPP
-#define COEFFICIENTS_HPP
-
-
-#include "mfem.hpp"
-#include "linalg/dual.hpp"
-
-namespace mfem {
-
-
-class StrainCoefficient:public VectorCoefficient{
-public:
-    StrainCoefficient(GridFunction& disp_):VectorCoefficient(disp_.VectorDim()*(disp_.VectorDim()+1)/2)
-    {
-        disp=&disp_;
-        g.SetSize(disp->VectorDim());
-        e.SetSize(disp->VectorDim());
-    }
-
-    virtual
-    ~StrainCoefficient(){}
-
-    virtual
-    void Eval(Vector &V, ElementTransformation &T,
-              const IntegrationPoint &ip)
-    {
-        T.SetIntPoint(&ip);
-        disp->GetVectorGradient(T,g);
-        for(int i=0;i<disp->VectorDim();i++){
-        for(int j=i;j<disp->VectorDim();j++){
-            e(i,j)=0.5*(g(i,j)+g(j,i));
-        }}
-
-        V.SetSize(GetVDim());
-
-        if(disp->VectorDim()==3){
-            V(0)=e(0,0);
-            V(1)=e(1,1);
-            V(2)=e(2,2);
-            V(3)=e(1,2);
-            V(4)=e(0,2);
-            V(5)=e(0,1);
-        }else{
-            V(0)=e(0,0);
-            V(1)=e(1,1);
-            V(2)=e(0,1);
-        }
-    }
-
-private:
-    GridFunction* disp;
-    DenseMatrix g;
-    DenseMatrix e;
-
-};
-
-
-class EngStrainCoefficient: public mfem::VectorCoefficient{
-public:
-
-    EngStrainCoefficient(GridFunction& disp_):VectorCoefficient(disp_.VectorDim()*(disp_.VectorDim()+1)/2)
-    {
-        disp=&disp_;
-        g.SetSize(disp->VectorDim());
-        e.SetSize(disp->VectorDim());
-    }
-
-    virtual
-    ~EngStrainCoefficient(){}
-
-    virtual
-    void Eval(Vector &V, ElementTransformation &T,
-              const IntegrationPoint &ip)
-    {
-        T.SetIntPoint(&ip);
-        disp->GetVectorGradient(T,g);
-        for(int i=0;i<disp->VectorDim();i++){
-        for(int j=i;j<disp->VectorDim();j++){
-            e(i,j)=0.5*(g(i,j)+g(j,i));
-        }}
-
-        V.SetSize(GetVDim());
-
-        if(disp->VectorDim()==3){
-            V(0)=e(0,0);
-            V(1)=e(1,1);
-            V(2)=e(2,2);
-            V(3)=2.0*e(1,2);
-            V(4)=2.0*e(0,2);
-            V(5)=2.0*e(0,1);
-        }else{
-            V(0)=e(0,0);
-            V(1)=e(1,1);
-            V(2)=2.0*e(0,1);
-        }
-    }
-
-
-private:
-    GridFunction* disp;
-    DenseMatrix g;
-    DenseMatrix e;
-
-
-};
-
-template<typename ElastMaterial>
-class StressCoefficient:public VectorCoefficient
-{
-public:
-    StressCoefficient(ElastMaterial& mat_, MatrixCoefficient& str_):
-        VectorCoefficient(str_.GetWidth()*(str_.GetWidth()+1)/2)
-    {
-        mat=&mat_;
-        str=&str_;
-        ss.SetSize(str_.GetWidth());
-
-    }
-
-    virtual
-    void Eval(Vector &V, ElementTransformation &T,
-              const IntegrationPoint &ip)
-    {
-        str->Eval(ss,T,ip);
-        mat->EvalStress(st,ss,T,ip);
-        if(str->GetWidth()==3){
-            V(0)=st(0,0);
-            V(1)=st(1,1);
-            V(2)=st(2,2);
-            V(3)=st(1,2);
-            V(4)=st(0,2);
-            V(5)=st(0,1);
-        }else{
-            V(0)=st(0,0);
-            V(1)=st(1,1);
-            V(2)=st(0,1);
-        }
-    }
-
-private:
-    ElastMaterial* mat;
-    MatrixCoefficient* str;
-
-    DenseMatrix ss; //strain
-    DenseMatrix st; //stress
-};
-
-template<typename ElastMaterial>
-class VonMisesStressCoefficient:public Coefficient
-{
-public:
-    VonMisesStressCoefficient(ElastMaterial& mat_, MatrixCoefficient& str_){
-        mat=&mat_;
-        str=&str_;//stress
-        ss.SetSize(str->GetWidth());
-        st.SetSize(str->GetWidth());
-    }
-
-    virtual
-    double Eval(ElementTransformation &T,
-                 const IntegrationPoint &ip)
-    {
-        str->Eval(ss,T,ip);
-        mat->EvalStress(st,ss,T,ip);
-
-        double res=0.0;
-
-        if(str->GetWidth()==3){
-            res=(st(0,0)-st(1,1))*(st(0,0)-st(1,1));
-            res=res+(st(1,1)-st(2,2))*(st(1,1)-st(2,2));
-            res=res+(st(0,0)-st(2,2))*(st(0,0)-st(2,2));
-            res=res+6*st(1,2)*st(1,2);
-            res=res+6*st(0,2)*st(0,2);
-            res=res+6*st(0,1)*st(0,1);
-        }else{
-            res=(st(0,0)-st(1,1))*(st(0,0)-st(1,1));
-            res=res+st(1,1)*st(1,1);
-            res=res+st(0,0)*st(0,0);
-            res=res+6*st(0,1)*st(0,1);
-        }
-
-        return sqrt(res);
-
-    }
-
-private:
-   ElastMaterial* mat;
-   MatrixCoefficient* str;
-   DenseMatrix st; //stress
-   DenseMatrix ss; //strain
-
-};
-
-class IsoElastMat
-{
-public:
-    IsoElastMat()
-    {
-        E=1.0;
-        nu=0.2;
-    }
-
-    void SetLameParam(double lam, double mu)
-    {
-        E=mu*(3.0*lam+2.0*mu)/(lam+mu);
-        nu=lam/(2.0*(lam+mu));
-    }
-
-    void SetElastParam(double E_,double nu_)
-    {
-        E=E_;
-        nu=nu_;
-    }
-
-    void SetE(double E_){ E=E_;}
-
-    void SetPoisson(double nu_){ nu=nu_;}
-
-
-    void EvalStress(DenseMatrix st, DenseMatrix ss,
-                    ElementTransformation &T,
-                    const IntegrationPoint &ip)
-    {
-        double cc=E/((1.0+nu)*(1.0-2.0*nu));
-        if(ss.Width()==3){
-            st(0,0)=cc*((1.0-nu)*ss(0,0)+nu*ss(1,1)+nu*ss(2,2));
-            st(1,1)=cc*(nu*ss(0,0)+(1.0-nu)*ss(1,1)+nu*ss(2,2));
-            st(2,2)=cc*(nu*ss(0,0)+nu*ss(1,1)+(1.0-nu)*ss(2,2));
-            st(1,2)=cc*(1.0-2.0*nu)*ss(1,2); st(2,1)=cc*(1.0-2.0*nu)*ss(2,1);
-            st(0,2)=cc*(1.0-2.0*nu)*ss(0,2); st(2,0)=cc*(1.0-2.0*nu)*ss(2,0);
-            st(0,1)=cc*(1.0-2.0*nu)*ss(0,1); st(1,0)=cc*(1.0-2.0*nu)*ss(1,0);
-        }else{
-            st(0,0)=cc*((1.0-nu)*ss(0,0)+nu*ss(1,1));
-            st(1,1)=cc*(nu*ss(0,0)+(1.0-nu)*ss(1,1));
-            //st(2,2)=cc*(nu*ss(0,0)+nu*ss(1,1));
-            st(0,1)=cc*(1.0-2.0*nu)*ss(0,1); st(1,0)=cc*(1.0-2.0*nu)*ss(1,0);
-        }
-    }
-
-    template<typename dtype>
-    void EvalStress(dtype* st, dtype* ss,
-                    ElementTransformation &T,
-                    const IntegrationPoint &ip)
-    {
-        double cc=E/((1.0+nu)*(1.0-2.0*nu));
-        st[0+3*0]=cc*((1.0-nu)*ss[0+3*0]+nu*ss[1+3*1]+nu*ss[2+3*2]);
-        st[1+3*1]=cc*(nu*ss[0+3*0]+(1.0-nu)*ss[1+3*1]+nu*ss[2+3*2]);
-        st[2+3*2]=cc*(nu*ss[0+3*0]+nu*ss[1+3*1]+(1.0-nu)*ss[2+3*2]);
-        st[1+3*2]=cc*(1.0-2.0*nu)*ss[1+3*2]; st[2+3*1]=cc*(1.0-2.0*nu)*ss[2+3*1];
-        st[0+3*2]=cc*(1.0-2.0*nu)*ss[0+3*2]; st[2+3*0]=cc*(1.0-2.0*nu)*ss[2+3*0];
-        st[0+3*1]=cc*(1.0-2.0*nu)*ss[0+3*1]; st[1+3*0]=cc*(1.0-2.0*nu)*ss[1+3*0];
-
-    }
-
-    // mat[9x9] ss[9]
-    template<typename dtype>
-    void EvalGrad(dtype* mat, dtype* ss,
-                  ElementTransformation &T,
-                  const IntegrationPoint &ip)
-    {
-        double cc=E/((1.0+nu)*(1.0-2.0*nu));
-        for(int i=0;i<81;i++){ mat[i]=dtype(0.0);}
-
-        mat[0*9 + 0]=cc*(1.0-nu);
-        mat[0*9 + 4]=cc*nu;
-        mat[0*9 + 8]=cc*nu;
-
-        mat[4*9 + 0]=cc*nu;
-        mat[4*9 + 4]=cc*(1.0-nu);
-        mat[4*9 + 8]=cc*nu;
-
-        mat[8*9+  0]=cc*nu;
-        mat[8*9+  4]=cc*nu;
-        mat[8*9+  8]=cc*(1.0-nu);
-
-        mat[1+3*2+(1+3*2)*9]=cc*(1.0-2.0*nu);
-        mat[0+3*2+(0+3*2)*9]=cc*(1.0-2.0*nu);
-        mat[0+3*1+(0+3*1)*9]=cc*(1.0-2.0*nu);
-
-        mat[2+3*1+(2+3*1)*9]=cc*(1.0-2.0*nu);
-        mat[2+3*0+(2+3*0)*9]=cc*(1.0-2.0*nu);
-        mat[1+3*0+(1+3*0)*9]=cc*(1.0-2.0*nu);
-    }
-
-
-    void EvalGrad(DenseMatrix& mat,
-                  Vector& ee, // strains
-                  ElementTransformation &T,
-                  const IntegrationPoint &ip)
-    {
-        mat.SetSize(9);
-        EvalGrad(mat.GetData(), ee.GetData(), T, ip);
-    }
-
-private:
-    double E;
-    double nu;
-};
-
-
-class J2YieldFunction
-{
-public:
-
-    J2YieldFunction()
-    {
-        sigma_0=1.0;
-        H=0.0;
-        beta=0.0;
-    }
-
-    J2YieldFunction(double ssy_,double H_, double beta_)
-    {
-        sigma_0=ssy_;
-        H=H_;
-        beta=beta_;
-    }
-
-    void Set(double ssy_,double H_, double beta_)
-    {
-        sigma_0=ssy_;
-        H=H_;
-        beta=beta_;
-    }
-
-    void SetJ2Param(double ssy_, double H_, double beta_)
-    {
-        sigma_0=ssy_;
-        H=H_;
-        beta=beta_;
-    }
-
-    //str[9] - stress tensor
-    //ip[2]:
-    //      ip[0] accumulated plastic strain (possitive)
-    //      ip[1] filtered accumulated plastic strain (possitive)
-    template<typename dtype0,typename dtype1,typename dtype2 >
-    dtype0 EvalI(dtype1* str,dtype2* ip)
-    {
-        dtype0 p=(str[0]+str[4]+str[8])/3.0;
-        dtype0 dstr[9];
-        for(int i=0;i<9;i++){ dstr[i]=str[i];}
-        dstr[0]=dstr[0]-p;
-        dstr[4]=dstr[4]-p;
-        dstr[8]=dstr[8]-p;
-        dtype0 se=dstr[0]*dstr[0];
-        for(int i=1;i<9;i++){ se=se+dstr[i]*dstr[i];}
-        se=sqrt(3.0*se/2.0);
-
-        dtype0 sy=sigma_0+H*ip[0];
-        dtype0 tv=exp(-beta*ip[1]);
-        //return se-tv*tv*sy*sy;
-        return se-tv*sy;
-    }
-
-    //return gradients of the yield function with respect to str and ip
-    template<typename dtype>
-    dtype EvalG(dtype* gstr, dtype* gip, dtype* str, dtype* ip)
-    {
-        dtype t1,t7,t9,t11,t13,t17,t19,t21,t23,t25,t26,t27,t29,t30,t40;
-        t1 = str[0]*str[0];
-        t7 = str[1]*str[1];
-        t9 = str[2]*str[2];
-        t11 = str[3]*str[3];
-        t13 = str[4]*str[4];
-        t17 = str[5]*str[5];
-        t19 = str[6]*str[6];
-        t21 = str[7]*str[7];
-        t23 = str[8]*str[8];
-        t25 = -4.0*str[0]*str[4]-4.0*str[0]*str[8]-4.0*str[4]*str[8]+4.0*t1+6.0*
-                t11+4.0*t13+6.0*t17+6.0*t19+6.0*t21+4.0*t23+6.0*t7+6.0*t9;
-        t26 = sqrt(t25);
-        t27 = 1.0/t26;
-        t29 = 4.0*str[4];
-        t30 = 4.0*str[8];
-        t40 = 4.0*str[0];
-        gstr[0] = t27*(8.0*str[0]-t29-t30)/4.0;
-        gstr[1] = 3.0*t27*str[1];
-        gstr[2] = 3.0*t27*str[2];
-        gstr[3] = 3.0*t27*str[3];
-        gstr[4] = t27*(-t40+8.0*str[4]-t30)/4.0;
-        gstr[5] = 3.0*t27*str[5];
-        gstr[6] = 3.0*t27*str[6];
-        gstr[7] = 3.0*t27*str[7];
-        gstr[8] = t27*(-t40-t29+8.0*str[8])/4.0;
-
-        gip[0] = -exp(-beta*ip[1])*H;
-        gip[1] = beta*exp(-beta*ip[1])*(H*ip[0]+sigma_0);
-
-        dtype sy=sigma_0+H*ip[0];
-        dtype tv=exp(-beta*ip[1]);
-        return t26/2.0-tv*sy;
-    }
-
-    double Eval(const Vector& str,const Vector& ivar)
-    {
-        return EvalI<double,double,double>(str.GetData(),ivar.GetData());
-    }
-
-    double Eval(double* str, double* ivar)
-    {
-        return EvalI<double,double,double>(str,ivar);
-    }
-
-    double EvalFGrad(Vector& drstr, Vector& divar,
-                    Vector& str, Vector& ivar)
-    {
-        //drstr.SetSize(str.Size());
-        //divar.SetSize(ivar.Size());
-        typedef internal::dual<double, double> ADFloatType;
-        ADFloatType tst[str.Size()];
-        ADFloatType var[ivar.Size()];
-
-        for(int i=0;i<str.Size();i++){
-            tst[i].value=str[i];
-            tst[i].gradient=0.0;
-        }
-
-        for(int i=0;i<ivar.Size();i++){
-            var[i].value=ivar[i];
-            var[i].gradient=0.0;
-        }
-
-        ADFloatType rez;
-        for(int i=0;i<str.Size();i++){
-            tst[i].gradient=1.0;
-            rez=EvalI<ADFloatType,ADFloatType,ADFloatType>(tst,var);
-            drstr[i]=rez.gradient;
-            tst[i].gradient=0.0;
-        }
-
-        for(int i=0;i<ivar.Size();i++){
-            var[i].gradient=1.0;
-            rez=EvalI<ADFloatType,ADFloatType,ADFloatType>(tst,var);
-            divar[i]=rez.gradient;
-            var[i].gradient=0.0;
-        }
-
-        return rez.value;
-    }
-
-private:
-    double H;
-    double beta;
-    double sigma_0;
-
-};
-
-
-class DruckerPragerYieldFunction
-{
-public:
-    DruckerPragerYieldFunction(double ssy_,double alpha_)
-    {
-        sigma_0=ssy_;
-        alpha=alpha_;
-    }
-
-    template<typename dtype0, typename dtype1, typename dtype2>
-    dtype0 Eval(dtype1* str,dtype2* ip)
-    {
-        dtype0 p=(str[0]+str[4]+str[8])/3.0;
-        dtype0 dstr[9];
-        for(int i=0;i<9;i++){ dstr[i]=str[i];}
-        dstr[0]=dstr[0]-p;
-        dstr[4]=dstr[1]-p;
-        dstr[8]=dstr[8]-p;
-        dtype0 J2=dstr[0]*dstr[0];
-        for(int i=1;i<9;i++){ J2=J2+dstr[i]*dstr[i];}
-        J2=J2/2.0;
-        return sqrt(J2)-3.0*alpha*p-sigma_0/sqrt(3.0);
-    }
-
-    double operator()(Vector& str,Vector& ivar)
-    {
-        return Eval<double>(str.GetData(),ivar.GetData());
-    }
-
-private:
-
-    double sigma_0;
-    double alpha;
-};
-
-class MatsuokaNakaiYieldFunction
-{
-public:
-    MatsuokaNakaiYieldFunction(double phi_)
-    {
-        phi=phi_;
-    }
-
-    template<typename dtype0, typename dtype1, typename dtype2>
-    dtype0 Eval(dtype1* str,dtype2* ip)
-    {
-        dtype0 I1=(str[0]+str[4]+str[8]);
-        dtype0 I2=str[0]*str[4]+str[4]*str[8]+str[8]*str[0]
-                              -str[3]*str[3]-str[6]*str[6]-str[7]*str[7];
-        dtype0 I3= str[0]*str[4]*str[8]+str[1]*str[5]*str[6]+str[3]*str[7]*str[2]
-                 -str[2]*str[4]*str[6]-str[1]*str[3]*str[8]-str[5]*str[7]*str[0];
-
-        double tp=tan(phi);
-        return I1*I2 - (9.0+8.0*tp*tp)*I3;
-    }
-
-    double operator()(Vector& str,Vector& ivar)
-    {
-        return Eval<double>(str.GetData(),ivar.GetData());
-    }
-private:
-    double phi;
-};
-
-
-
-
-
-template<typename elmat, typename yfunc>
-class StressEval
-{
-public:
-
-    StressEval(elmat* mat_, yfunc* yf_, double lerr_=1e-8){
-        mat=mat_;
-        yf=yf_;
-        lerr=lerr_;
-    }
-
-    void SetPlasticStrain(const Vector& ep_)
-    {
-        for(int i=0;i<9;i++){
-            epn[i]=ep_[i];
-        }
-    }
-
-    void SetInternalParameters(const Vector& ip_)
-    {
-        ipn[0]=ip_[0];
-        ipn[1]=ip_[1];
-    }
-
-    void SetStrain(const Vector &ee_)
-    {
-        for(int i=0;i<9;i++){
-            ee[i]=ee_[i];
-        }
-    }
-
-    void EvalTangent(DenseMatrix& Cep,
-                     ElementTransformation &T,
-                     const IntegrationPoint &ip)
-    {
-        Vector ees(9);
-        Vector css(9);
-        for(int i=0;i<9;i++){
-            ees[i]=ee[i]-epn[i];
-        }
-
-        mat->EvalStress(css.GetData(),ees.GetData(),T,ip);
-        mat->EvalGrad(Cep,ees,T,ip);
-        double f=yf->Eval(css.GetData(),ipn);
-        if(f<0.0){
-            return;}
-        Vector dip(2);
-        //elasto-plastic behaviour - modify the elastic tensor
-        Vector r(9);
-        //yf->EvalFGrad(r,dip,css,ipn);
-        yf->EvalG(r.GetData(), dip.GetData(), css.GetData(),ipn);
-        double H;
-
-        mat->EvalStress(css.GetData(),r.GetData(),T,ip);
-        H=(css*r)-dip[0]*sqrt(2.0*(r*r)/3.0);
-        for(int i=0;i<9;i++){
-        for(int j=0;j<9;j++){
-            Cep(i,j)=Cep(i,j)-css(i)*css(j)/H;
-        }}
-
-    }
-
-    /// vin[stress[9], strainp[9], ip[2], multiplier[1]]
-    template<typename dtype>
-    void EvalResidual(dtype* rr, dtype* vin,
-                      ElementTransformation &T,
-                      const IntegrationPoint &ip)
-    {
-        dtype* css=vin;   //current stress
-        dtype* cep=vin+9; //current plastic strain
-        dtype* cip=vin+18;//current internal variables
-        dtype* lam=vin+20;//current Largrange multiplier
-
-        dtype* rr1=rr;
-        dtype* rr2=rr+9;
-        dtype* rr3=rr+18;
-        dtype* rr4=rr+20;
-
-        dtype tv[9];
-        dtype hh[2];
-        dtype vv;
-
-        //stress residual
-        for(int i=0;i<9;i++){
-            tv[i]=ee[i]-cep[i];
-        }
-        mat->EvalStress(rr1,tv,T,ip);
-        for(int i=0;i<9;i++){
-            rr1[i]=css[i]-rr1[i];
-        }
-
-        rr4[0]=yf->EvalG(tv, hh, css, cip);
-        for(int i=0;i<9;i++){
-            rr2[i]=cep[i]-epn[i]-(*lam)*tv[i];
-        }
-
-        //internal parameters
-        vv=tv[0]*tv[0];
-        for(int i=1;i<9;i++){ vv=vv+tv[i]*tv[i];}
-        rr3[0]=cip[0]-ipn[0]-(*lam)*sqrt((2.0/3.0)*vv);
-        rr3[1]=cip[1]-ipn[1];
-    }
-
-    //E
-    void EvalNewton(DenseMatrix& tmat, Vector& rr,
-                    Vector& vv,
-                    ElementTransformation &T,
-                    const IntegrationPoint &ip){
-        typedef internal::dual<double, double> ADFloatType;
-        ADFloatType tv[21];
-        ADFloatType cv[21];
-        for(int i=0;i<21;i++){
-            cv[i].value=vv[i];
-            cv[i].gradient=0.0;
-        }
-
-        for(int i=0;i<21;i++){
-            cv[i].gradient=1.0;
-            EvalResidual(tv, cv,T,ip);
-            for(int j=0;j<21;j++){
-                tmat(j,i)=tv[j].gradient;
-            }
-            cv[i].gradient=0.0;
-        }
-
-        for(int i=0;i<21;i++){
-            rr[i]=tv[i].value;
-        }
-    }
-
-    void Solve(Vector& css, Vector& cep, Vector& cip,
-               ElementTransformation &T,
-               const IntegrationPoint &ip)
-    {
-        Vector ees(9);
-        for(int i=0;i<9;i++){
-            ees[i]=ee[i]-epn[i];
-            cep[i]=epn[i];
-        }
-
-        cip[0]=ipn[0];
-        cip[1]=ipn[1];
-
-        mat->EvalStress(css.GetData(),ees.GetData(),T,ip);
-
-        double f=yf->Eval(css,cip);
-        if(f<0.0){return;}
-        SolveEP(css,cep,cip,T,ip);
-    }
-
-    void SolveEP(Vector& css, Vector& cep, Vector& cip,
-               ElementTransformation &T,
-               const IntegrationPoint &ip, double err=1e-8){
-        Vector vv(21);
-        Vector rr(21);
-        Vector tv(21);
-        DenseMatrix mm(21);
-        for(int i=0;i<9;i++){
-            vv[i]=css[i];
-            vv[i+9]=cep[i];
-        }
-        vv[18]=cip[0];
-        vv[19]=cip[1];
-        vv[20]=0.0;
-
-        EvalNewton(mm,rr,vv,T,ip);
-
-        /*
-        std::cout<<std::endl;
-        mm.PrintMatlab(std::cout);
-        std::cout<<"rr="<<std::endl;
-        rr.Print(std::cout,21);
-        */
-
-        double cerr=rr.Norml2();
-        int it=0;
-        while(cerr>err){
-            DenseMatrixInverse im(mm);
-            im.Mult(rr,tv);
-            vv.Add(-1.0,tv);
-            EvalNewton(mm,rr,vv,T,ip);
-            cerr=rr.Norml2();
-            //std::cout<<"it:"<<it<<" err="<<cerr<<std::endl;
-            it++;
-            if(it>100){
-                mfem_error("Maximum number (100) of iterations have been"
-                           "reached in the elasto-plastic material update! \n");
-                break;}
-        }
-
-        for(int i=0;i<9;i++){
-            css[i]=vv[i];
-            cep[i]=vv[i+9];
-        }
-        cip[0]=vv[18];
-        cip[1]=vv[19];
-    }
-
-
-
-private:
-    elmat* mat;
-    yfunc* yf;
-
-    double ee[9]; //total strain
-    double epn[9]; //plastic strain
-    double ipn[2]; //internal variables
-
-    double lerr;
-
-};
-
-
-}
-
-#endif

miniapps/elplast/elpl_test.cpp

@@ -1,331 +0,0 @@
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-#include "coefficients.hpp"
-#include "elplast.hpp"
-
-using namespace std;
-using namespace mfem;
-
-
-class MyVectorCoeff:public VectorCoefficient
-{
-public:
-    MyVectorCoeff():VectorCoefficient(3)
-    {
-        SetTime(1.0);
-    }
-
-    virtual
-    ~MyVectorCoeff(){}
-
-
-    void Eval(Vector &V, ElementTransformation &T,
-               const IntegrationPoint &ip)
-    {
-        Vector xx(3);
-        T.Transform(ip,xx);
-        V=xx; V*=0.0;
-        V[0]=xx[0]*GetTime();
-        //V[1]=xx[0]*GetTime();
-        //V[2]=xx[2]*GetTime();
-    }
-
-private:
-};
-
-class MyForceCoeff:public VectorCoefficient
-{
-public:
-    MyForceCoeff():VectorCoefficient(3)
-    {
-        SetTime(1.0);
-    }
-
-    virtual
-    ~MyForceCoeff(){}
-
-
-    void Eval(Vector &V, ElementTransformation &T,
-               const IntegrationPoint &ip)
-    {
-        Vector xx(3);
-        T.Transform(ip,xx);
-        V=xx; V*=0.0;
-        if(xx[0]>1.9){V[0]=GetTime();}
-        //V[1]=xx[1]*GetTime();
-        //V[2]=xx[2]*GetTime();
-    }
-
-private:
-};
-
-
-class MyBlockPrec:public BlockDiagonalPreconditioner
-{
-public:
-    MyBlockPrec(const Array<int> & offsets, ParFiniteElementSpace* fesel):BlockDiagonalPreconditioner(offsets)
-    {
-        pr1=new HypreBoomerAMG();
-        pr1->SetPrintLevel(0);
-        pr1->SetElasticityOptions(fesel);
-        pr2=new HypreBoomerAMG();
-        pr2->SetPrintLevel(0);
-
-    }
-
-    virtual void SetOperator(const Operator &op)
-    {
-        const BlockOperator& blo=static_cast<const BlockOperator&>(op);
-        pr1->SetOperator(blo.GetBlock(0,0));
-        pr2->SetOperator(blo.GetBlock(1,1));
-        SetDiagonalBlock(0,pr1);
-        SetDiagonalBlock(1,pr2);
-    }
-
-    virtual ~MyBlockPrec()
-    {
-        delete pr1;
-        delete pr2;
-    }
-
-private:
-    HypreBoomerAMG* pr1;
-    HypreBoomerAMG* pr2;
-
-};
-
-
-
-
-int main(int argc, char *argv[])
-{
-    // 0. Initialize MPI and HYPRE.
-    Mpi::Init(argc, argv);
-    int num_procs = Mpi::WorldSize();
-    int myid = Mpi::WorldRank();
-    Hypre::Init();
-
-    // 1. Parse command-line options.
-    const char *mesh_file = "../data/beam-tri.mesh";
-    int serial_ref_levels = 0;
-    int order = 2;
-    //bool static_cond = false;
-    bool static_cond = true;
-    bool visualization = 1;
-
-    OptionsParser args(argc, argv);
-    args.AddOption(&mesh_file, "-m", "--mesh",
-                   "Mesh file to use.");
-    args.AddOption(&serial_ref_levels, "-rs", "--refine-serial",
-                   "Number of uniform serial refinements (before parallel"
-                   " partitioning)");
-    args.AddOption(&order, "-o", "--order",
-                   "Finite element order (polynomial degree).");
-    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                   "--no-static-condensation", "Enable static condensation.");
-    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                   "--no-visualization",
-                   "Enable or disable GLVis visualization.");
-    args.Parse();
-    if (!args.Good())
-    {
-       if (myid == 0)
-       {
-          args.PrintUsage(cout);
-       }
-       return 1;
-    }
-    if (myid == 0)
-    {
-       args.PrintOptions(cout);
-    }
-
-    // 2. Read the mesh from the given mesh file. We can handle triangular,
-    //    quadrilateral, tetrahedral, and hexahedral meshes with the same code.
-    Mesh mesh(mesh_file, 1, 1);
-    int dim = mesh.Dimension();
-    MFEM_VERIFY(mesh.SpaceDimension() == dim, "invalid mesh");
-
-    // 3. Refine the mesh before parallel partitioning. Since a NURBS mesh can
-    //    currently only be refined uniformly, we need to convert it to a
-    //    piecewise-polynomial curved mesh. First we refine the NURBS mesh a bit
-    //    more and then project the curvature to quadratic Nodes.
-    if (mesh.NURBSext && serial_ref_levels == 0)
-    {
-       serial_ref_levels = 2;
-    }
-    for (int i = 0; i < serial_ref_levels; i++)
-    {
-       mesh.UniformRefinement();
-    }
-    if (mesh.NURBSext)
-    {
-       mesh.SetCurvature(2);
-    }
-    mesh.EnsureNCMesh();
-
-    ParMesh pmesh(MPI_COMM_WORLD, mesh);
-    mesh.Clear();
-
-    // 4. Define a finite element space on the mesh. The polynomial order is
-    //    one (linear) by default, but this can be changed on the command line.
-    H1_FECollection fec(order, dim);
-    ParFiniteElementSpace ufespace(&pmesh, &fec, dim);
-    ParFiniteElementSpace efespace(&pmesh, &fec, 1);
-
-    QuadratureSpace  qfes(&pmesh,3*order);
-    QuadratureFunction kappa;
-    QuadratureFunction eep;
-
-    kappa.SetSpace(&qfes,1); kappa=0.0;
-    eep.SetSpace(&qfes,6);   eep=0.0;
-
-    ParGridFunction u(&ufespace); u=0.0;
-    ParGridFunction ep(&efespace); ep=0.0;
-    MyVectorCoeff co; co.SetTime(0.5);
-    MyForceCoeff  fo; fo.SetTime(0.5);
-
-    u.ProjectCoefficient(co);
-
-    Array<ParFiniteElementSpace *> pfes(2);
-    pfes[0]=&ufespace;
-    pfes[1]=&efespace;
-    ParBlockNonlinearForm* bnl=new ParBlockNonlinearForm();
-    bnl->SetParSpaces(pfes);
-
-    //add the integrator
-    NLElPlastIntegrator* itgr=new NLElPlastIntegrator(1,0.2,1.0,0.01);
-    itgr->SetPlasticStrains(eep,kappa);
-    //itgr->SetForce(fo);
-    bnl->AddDomainIntegrator(itgr);
-
-    Array<int> offset;
-    BlockVector bv; bv.Update(bnl->GetBlockTrueOffsets()); bv=0.0;
-    BlockVector rv; rv.Update(bnl->GetBlockTrueOffsets()); rv=0.0;
-
-    u.GetTrueDofs(bv.GetBlock(0));
-
-    bnl->Mult(bv,rv);
-
-    BlockOperator& op=bnl->GetGradient(bv);
-
-
-    // Define the essential boundary attributes
-    Array<int> ess_bdr_u(pmesh.bdr_attributes.Max());
-    Array<int> ess_bdr_e(pmesh.bdr_attributes.Max());
-    Array<Array<int>*> ess_bdr(2);
-    ess_bdr[0]=&ess_bdr_u;
-    ess_bdr[1]=&ess_bdr_e;
-
-    ess_bdr_u = 1; ess_bdr_u[0] = 1;
-    ess_bdr_e = 0;
-
-    Array<Vector*> nrhs(2);
-    nrhs[0]=nullptr;
-    nrhs[1]=nullptr;
-
-    bnl->SetEssentialBC(ess_bdr,nrhs);
-
-
-    ParGridFunction du(u);
-    du.SetFromTrueDofs(bv.GetBlock(0));
-    ep.SetFromTrueDofs(bv.GetBlock(1));
-
-    // ParaView output.
-    ParaViewDataCollection dacol("ParaView", &pmesh);
-    dacol.SetLevelsOfDetail(order);
-    dacol.SetDataFormat(VTKFormat::ASCII);
-    dacol.RegisterField("odisp", &u);
-    dacol.RegisterField("ddisp", &du);
-    dacol.RegisterField("ep",&ep);
-
-    /*
-    dacol.SetTime(1.0);
-    dacol.SetCycle(1);
-    dacol.Save();
-    */
-
-
-
-
-    /*
-    BlockDiagonalPreconditioner* blprec=new BlockDiagonalPreconditioner(bnl->GetBlockTrueOffsets());
-    HypreBoomerAMG* pr1=new HypreBoomerAMG();
-    pr1->SetPrintLevel(1);
-    HypreBoomerAMG* pr2=new HypreBoomerAMG();
-    pr2->SetPrintLevel(1);
-    */
-
-
-    MyBlockPrec* blprec=new MyBlockPrec(bnl->GetBlockTrueOffsets(), &ufespace);
-
-    /*
-    GMRESSolver* gmres = new GMRESSolver(MPI_COMM_WORLD);
-    gmres->SetAbsTol(1e-12);
-    gmres->SetRelTol(1e-9);
-    gmres->SetMaxIter(500);
-    gmres->SetPrintLevel(2);
-    gmres->SetPreconditioner(*blprec);
-    */
-
-    CGSolver* cg = new CGSolver(MPI_COMM_WORLD);
-    cg->SetAbsTol(1e-12);
-    cg->SetRelTol(1e-9);
-    cg->SetMaxIter(500);
-    cg->SetPrintLevel(2);
-    cg->SetPreconditioner(*blprec);
-
-
-    NewtonSolver* ns=new  NewtonSolver(MPI_COMM_WORLD);
-
-    ns->iterative_mode = true;
-    ns->SetSolver(*cg);
-    ns->SetOperator(*bnl);
-    ns->SetPrintLevel(1);
-    ns->SetRelTol(1e-5);
-    ns->SetAbsTol(1e-12);
-    ns->SetMaxIter(8);
-
-    u.GetTrueDofs(bv.GetBlock(0));
-    //bv=0.0;
-    Vector b;
-    //ns->Mult(b,bv);
-
-    for(int bi=0;bi<10;bi++){
-        co.SetTime(0.4+bi*0.1);
-        u.ProjectCoefficient(co);
-        u.GetTrueDofs(bv.GetBlock(0));
-        ns->Mult(b,bv);
-
-        //update the plasticity vars
-        itgr->SetUpdateFlag(true);
-        bnl->Mult(bv,rv);
-        itgr->SetUpdateFlag(false);
-
-        u.SetFromTrueDofs(bv.GetBlock(0));
-        ep.SetFromTrueDofs(bv.GetBlock(1));
-
-        dacol.SetTime(1.0+bi);
-        dacol.SetCycle(1+bi);
-        dacol.Save();
-
-    }
-
-
-    delete ns;
-    delete cg;
-    delete blprec;
-
-    MPI_Barrier(MPI_COMM_WORLD);
-
-
-
-    delete bnl;
-
-
-
-
-    return 0;
-}

miniapps/elplast/elplast.cpp

@@ -1,904 +0,0 @@
-#include "elplast.hpp"
-
-#include <fstream>
-#include <iostream>
-
-
-namespace mfem{
-
-
-ElPlastSolver::ElPlastSolver(mfem::ParMesh* mesh_,int vorder, int forder)
-{
-
-    pmesh=mesh_;
-
-    int dim=pmesh->Dimension();
-    vfec=new H1_FECollection(vorder,dim);
-    ffec=new H1_FECollection(forder,dim);
-
-    vfes=new ParFiniteElementSpace(pmesh,vfec,dim,Ordering::byVDIM);
-    ffes=new ParFiniteElementSpace(pmesh,ffec);
-
-    fdisp.SetSpace(vfes); fdisp=0.0;
-    adisp.SetSpace(vfes); adisp=0.0;
-
-
-    qfes=new QuadratureSpace(pmesh,3*std::max(vorder,forder));
-    kappa.SetSpace(qfes,1); kappa=0.0;
-    eep.SetSpace(qfes,6);   eep=0.0;
-    eee.SetSpace(qfes,6); eee=0.0;
-
-    nf=nullptr;
-    Array<mfem::ParFiniteElementSpace*> pf;
-    pf.Append(vfes);
-    pf.Append(ffes);
-
-    nf=new ParBlockNonlinearForm(pf);
-    rhs.Update(nf->GetBlockTrueOffsets()); rhs=0.0;
-    sol.Update(nf->GetBlockTrueOffsets()); sol=0.0;
-    adj.Update(nf->GetBlockTrueOffsets()); adj=0.0;
-
-
-    SetNewtonSolver();
-    SetLinearSolver();
-
-}
-
-
-ElPlastSolver::~ElPlastSolver()
-{
-
-    delete nf;
-
-    delete vfes;
-    delete ffes;
-    delete ffec;
-    delete vfec;
-
-    //delete the local forces
-    for(auto it=lvforce.begin();it!=lvforce.end();it++)
-    {
-        delete it->second;
-    }
-}
-
-void ElPlastSolver::SetNewtonSolver(double rtol, double atol,int miter, int prt_level)
-{
-    rel_tol=rtol;
-    abs_tol=atol;
-    max_iter=miter;
-    print_level=prt_level;
-}
-
-void ElPlastSolver::SetLinearSolver(double rtol, double atol, int miter)
-{
-    linear_rtol=rtol;
-    linear_atol=atol;
-    linear_iter=miter;
-}
-
-void ElPlastSolver::AddDispBC(int id, int dir, double val)
-{
-    if(dir==0){
-        bcx[id]=mfem::ConstantCoefficient(val);
-        AddDispBC(id,dir,bcx[id]);
-    }
-    if(dir==1){
-        bcy[id]=mfem::ConstantCoefficient(val);
-        AddDispBC(id,dir,bcy[id]);
-
-    }
-    if(dir==2){
-        bcz[id]=mfem::ConstantCoefficient(val);
-        AddDispBC(id,dir,bcz[id]);
-    }
-    if(dir==4){
-        bcx[id]=mfem::ConstantCoefficient(val);
-        bcy[id]=mfem::ConstantCoefficient(val);
-        bcz[id]=mfem::ConstantCoefficient(val);
-        AddDispBC(id,0,bcx[id]);
-        AddDispBC(id,1,bcy[id]);
-        AddDispBC(id,2,bcz[id]);
-    }
-}
-
-void ElPlastSolver::AddDispBC(int id, int dir, Coefficient &val)
-{
-    if(dir==0){ bccx[id]=&val; }
-    if(dir==1){ bccy[id]=&val; }
-    if(dir==2){ bccz[id]=&val; }
-    if(dir==4){ bccx[id]=&val; bccy[id]=&val; bccz[id]=&val;}
-    if(pmesh->Dimension()==2)
-    {
-        bccz.clear();
-    }
-}
-
-void ElPlastSolver::AddDispBC(int id, mfem::VectorCoefficient& val)
-{
-    bcca[id]=&val;
-}
-
-void ElPlastSolver::AddVolForce(int id, double fx, double fy, double fz)
-{
-    //check if the force id is already in the database
-    if(lvforce.find(id)!=lvforce.end()){
-        delete lvforce[id];
-    }
-    Vector tv(3); tv[0]=fx; tv[1]=fy; tv[2]=fz;
-    lvforce[id]=new mfem::VectorConstantCoefficient(tv);
-    volforce[id]=lvforce[id];
-}
-
-/// Adds vol force
-void ElPlastSolver::AddVolForce(int id, mfem::VectorCoefficient& ff)
-{
-    volforce[id]=&ff;
-}
-
-void ElPlastSolver::FSolve()
-{
-    // Set the BC
-    ess_tdofv.DeleteAll();
-    Array<int> ess_tdofx;
-    Array<int> ess_tdofy;
-    Array<int> ess_tdofz;
-
-    int dim=pmesh->Dimension();
-    {
-        for(auto it=bccx.begin();it!=bccx.end();it++)
-        {
-            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
-            ess_bdr=0;
-            ess_bdr[it->first -1]=1;
-            mfem::Array<int> ess_tdof_list;
-            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,0);
-            ess_tdofx.Append(ess_tdof_list);
-
-            mfem::VectorArrayCoefficient pcoeff(dim);
-            pcoeff.Set(0, it->second, false);
-            fdisp.ProjectBdrCoefficient(pcoeff, ess_bdr);
-        }
-
-        //copy tdofsx from velocity grid function
-        {
-            fdisp.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
-            for(int ii=0;ii<ess_tdofx.Size();ii++)
-            {
-                sol.GetBlock(0)[ess_tdofx[ii]]=rhs.GetBlock(0)[ess_tdofx[ii]];
-            }
-        }
-        ess_tdofv.Append(ess_tdofx);
-
-        for(auto it=bccy.begin();it!=bccy.end();it++)
-        {
-            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
-            ess_bdr=0;
-            ess_bdr[it->first -1]=1;
-            mfem::Array<int> ess_tdof_list;
-            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,1);
-            ess_tdofy.Append(ess_tdof_list);
-
-            mfem::VectorArrayCoefficient pcoeff(dim);
-            pcoeff.Set(1, it->second, false);
-            fdisp.ProjectBdrCoefficient(pcoeff, ess_bdr);
-        }
-        //copy tdofsy from velocity grid function
-        {
-            fdisp.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
-            for(int ii=0;ii<ess_tdofy.Size();ii++)
-            {
-                sol.GetBlock(0)[ess_tdofy[ii]]=rhs.GetBlock(0)[ess_tdofy[ii]];
-            }
-        }
-        ess_tdofv.Append(ess_tdofy);
-
-        if(dim==3){
-            for(auto it=bccz.begin();it!=bccz.end();it++)
-            {
-                mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
-                ess_bdr=0;
-                ess_bdr[it->first -1]=1;
-                mfem::Array<int> ess_tdof_list;
-                vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list,2);
-                ess_tdofz.Append(ess_tdof_list);
-
-                mfem::VectorArrayCoefficient pcoeff(dim);
-                pcoeff.Set(2, it->second, false);
-                fdisp.ProjectBdrCoefficient(pcoeff, ess_bdr);
-            }
-
-            //copy tdofsz from velocity grid function
-            {
-                fdisp.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
-                for(int ii=0;ii<ess_tdofz.Size();ii++)
-                {
-                    sol[ess_tdofz[ii]]=rhs[ess_tdofz[ii]];
-                }
-            }
-            ess_tdofv.Append(ess_tdofz);
-        }
-
-        //set vector coefficients
-        for(auto it=bcca.begin();it!=bcca.end();it++)
-        {
-            mfem::Array<int> ess_bdr(pmesh->bdr_attributes.Max());
-            ess_bdr=0;
-            ess_bdr[it->first -1]=1;
-            mfem::Array<int> ess_tdof_list;
-            vfes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
-            fdisp.ProjectBdrCoefficient(*(it->second), ess_bdr);
-            //copy tdofs from velocity grid function
-            fdisp.GetTrueDofs(rhs.GetBlock(0)); // use the rhs vector as a tmp vector
-            for(int ii=0;ii<ess_tdof_list.Size();ii++)
-            {
-                sol[ess_tdof_list[ii]]=rhs[ess_tdof_list[ii]];
-            }
-            ess_tdofv.Append(ess_tdof_list);
-        }
-    }
-
-
-
-}
-
-
-
-void NLElPlastIntegrator::AssembleElementVector(const Array<const FiniteElement *> &el,
-                                                ElementTransformation &Tr,
-                                                const Array<const Vector *> &elfun,
-                                                const Array<Vector *> &elvec)
-{
-    //the integrator works only for 3 dimensional problems
-    int dof_u = el[0]->GetDof();
-    int dof_e = el[1]->GetDof();
-
-    int dim =  Tr.GetSpaceDim();
-
-    elvec[0]->SetSize(dim*dof_u);
-    elvec[1]->SetSize(dof_e);
-
-    if (dim != 3)
-    {
-        mfem::mfem_error("NLElPlastIntegrator::AssembleElementVector"
-                                " is not defined on manifold meshes");
-    }
-
-    Vector uu(elfun[0]->GetData()+0*dof_u, dof_u);
-    Vector vv(elfun[0]->GetData()+1*dof_u, dof_u);
-    Vector ww(elfun[0]->GetData()+2*dof_u, dof_u);
-
-    Vector ru(elvec[0]->GetData()+0*dof_u, dof_u); ru=0.0;
-    Vector rv(elvec[0]->GetData()+1*dof_u, dof_u); rv=0.0;
-    Vector rw(elvec[0]->GetData()+2*dof_u, dof_u); rw=0.0;
-
-    Vector ep(elfun[1]->GetData(), dof_e);
-    Vector rp(elvec[1]->GetData(), dof_e); rp=0.0;
-
-    // temp storages for vectors and matrices
-    Vector su(dof_u); //shape functions for displacements
-    DenseMatrix du(dof_u,dim); //gradients of the shape functions
-    Vector dux; dux.SetDataAndSize(du.GetData()+0*dof_u,dof_u);
-    Vector duy; duy.SetDataAndSize(du.GetData()+1*dof_u,dof_u);
-    Vector duz; duz.SetDataAndSize(du.GetData()+2*dof_u,dof_u);
-    Vector se(dof_e); //shape functions for plastic strains
-    DenseMatrix de(dof_e,dim); //gradients of the shape functions
-    Vector grade(dim);
-
-    DenseMatrix vpss; //plastic strain at the integration points
-    DenseMatrix vkap; //accumulated plastic strain at the integration points
-    eep->GetElementValues(Tr.ElementNo,vpss);
-    kappa->GetElementValues(Tr.ElementNo,vkap);
-
-
-    DenseMatrix estrain(3,3);
-    Vector vestrain; vestrain.SetDataAndSize(estrain.GetData(),9);
-    DenseMatrix pstrain(3,3);
-    Vector vpstrain; vpstrain.SetDataAndSize(pstrain.GetData(),9);
-    DenseMatrix gradu(3,3); gradu=0.0;
-    Vector tv;
-
-    const IntegrationRule& ir=eep->GetElementIntRule(Tr.ElementNo);
-
-    double H=0.001;
-    double beta=0.0;
-    Vector vip(2); vip=0.0;
-
-    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> seval(&mat,&yf);
-    Vector stress(9);stress=0.0;
-
-    double w;
-    for (int i = 0; i < ir.GetNPoints(); i++)
-    {
-        const IntegrationPoint &ip = ir.IntPoint(i);
-        Tr.SetIntPoint(&ip);
-        w = Tr.Weight();
-        w = ip.weight * w;
-
-
-        //calculate the total strain
-        el[0]->CalcPhysShape(Tr,su);
-        el[0]->CalcPhysDShape(Tr,du);
-        tv.SetDataAndSize(gradu.GetData()+0*3,dim);
-        du.MultTranspose(uu,tv);
-        tv.SetDataAndSize(gradu.GetData()+1*3,dim);
-        du.MultTranspose(vv,tv);
-        tv.SetDataAndSize(gradu.GetData()+2*3,dim);
-        du.MultTranspose(ww,tv);
-        //compute the strain tensor
-        for(int ii=0;ii<dim;ii++){
-            for(int jj=ii+1;jj<dim;jj++){
-                estrain(ii,jj)=0.5*(gradu(ii,jj)+gradu(jj,ii));
-                estrain(jj,ii)=estrain(ii,jj);
-            }
-            estrain(ii,ii)=gradu(ii,ii);
-        }
-
-        //estrain.PrintMatlab(std::cout);
-
-        //set current plastic strain Voight indexing
-        {
-            pstrain(0,0)=vpss(0,i);
-            pstrain(1,1)=vpss(1,i);
-            pstrain(2,2)=vpss(2,i);
-            pstrain(1,2)=pstrain(2,1)=vpss(3,i);
-            pstrain(0,2)=pstrain(2,0)=vpss(4,i);
-            pstrain(0,1)=pstrain(1,0)=vpss(5,i);
-        }
-
-        //calculate the filtered plastic strain
-        el[1]->CalcPhysShape(Tr,se);
-        el[1]->CalcPhysDShape(Tr,de);
-        double epf=se*ep;//filtered accumulated plastic strain
-        double rry=ll->Eval(Tr,ip); //filter radius
-
-        //evaluate the material behaviour at the integration point
-        {
-            double Em=E->Eval(Tr,ip);
-            double nm=nu->Eval(Tr,ip);
-            double ssy=ss_y->Eval(Tr,ip);
-            mat.SetE(Em);
-            mat.SetPoisson(nm);
-            yf.Set(ssy,H,beta);
-
-            vip[0]=vkap(0,i); //set the accumulated plastic starin
-            vip[1]=epf;   //set the filtered plastic strain
-            seval.SetStrain(vestrain);
-            seval.SetPlasticStrain(vpstrain);
-            seval.SetInternalParameters(vip);
-            seval.Solve(stress,vpstrain,vip,Tr,ip);
-
-            if(flag_update){
-                vpss(0,i)=pstrain(0,0);
-                vpss(1,i)=pstrain(1,1);
-                vpss(2,i)=pstrain(2,2);
-                vpss(3,i)=pstrain(1,2);
-                vpss(4,i)=pstrain(0,2);
-                vpss(5,i)=pstrain(0,1);
-
-                vkap(0,i)=vip[0];
-            }
-
-            //stress.Print(std::cout,3);
-            //std::cout<<std::endl;
-        }
-
-        //assemble the RHS for displacements
-        ru.Add(stress[0]*w,dux);
-        ru.Add(0.5*stress[1]*w,duy); rv.Add(0.5*stress[1]*w,dux);
-        ru.Add(0.5*stress[2]*w,duz); rw.Add(0.5*stress[2]*w,dux);
-        ru.Add(0.5*stress[3]*w,duy); rv.Add(0.5*stress[3]*w,dux);
-        rv.Add(stress[4]*w,duy);
-        rv.Add(0.5*stress[5]*w,duz); rw.Add(0.5*stress[5]*w,duy);
-        ru.Add(0.5*stress[6]*w,duz); rw.Add(0.5*stress[6]*w,dux);
-        rv.Add(0.5*stress[7]*w,duz); rw.Add(0.5*stress[7]*w,duy);
-        rw.Add(stress[8]*w,duz);
-
-
-        //assemble the RHS for the accumulated plastic strain
-        rp.Add(-vip(0)*w,se);
-        rp.Add(epf*w,se);
-        de.MultTranspose(ep,grade);
-        tv.SetDataAndSize(de.GetData()+0*dof_e,dof_e);
-        rp.Add(rry*rry*grade[0]*w,tv);
-        tv.SetDataAndSize(de.GetData()+1*dof_e,dof_e);
-        rp.Add(rry*rry*grade[1]*w,tv);
-        tv.SetDataAndSize(de.GetData()+2*dof_e,dof_e);
-        rp.Add(rry*rry*grade[2]*w,tv);
-    }
-
-    //add the force
-    /*
-    if(force)
-    {
-        Vector lv(3*dof_u);
-        VectorDomainLFIntegrator li(*force);
-        li.AssembleRHSElementVect(*el[0],Tr,lv);
-        elvec[0]->Add(-1.0,lv);
-    }*/
-
-    //check the linear elasticity
-
-    /*
-    {
-        const IntegrationPoint &ip = ir.IntPoint(0);
-        double Em=E->Eval(Tr,ip);
-        double num=nu->Eval(Tr,ip);
-        double la=Em*num/((1.0+num)*(1.0-2.0*num));
-        double mu=Em/(2.0*(1+num));
-        ConstantCoefficient lc(la);
-        ConstantCoefficient mc(mu);
-        ElasticityIntegrator eli(lc,mc);
-        DenseMatrix K(dim*dof_u);
-        eli.AssembleElementMatrix(*el[0],Tr,K);
-        Vector rr(dim*dof_u);
-        //K.Mult(*(elfun[0]),rr);
-        //rr.Add(-1.0,*(elvec[0]));
-        *(elvec[0])=0.0;
-        K.Mult(*(elfun[0]),*(elvec[0]));
-        //std::cout<<"|rr|="<<rr.Norml2()<<" "<<elvec[0]->Norml2()<<std::endl;
-    }
-    */
-
-
-
-}
-
-void NLElPlastIntegrator::AssembleElementGrad(const Array<const FiniteElement *> &el,
-                                              ElementTransformation &Tr,
-                                              const Array<const Vector *> &elfun,
-                                              const Array2D<DenseMatrix *> &elmat)
-{
-    //the integrator works only for 3 dimensional problems
-    int dof_u = el[0]->GetDof();
-    int dof_e = el[1]->GetDof();
-
-    elmat(0,0)->SetSize(3*dof_u,3*dof_u); (*elmat(0,0))=0.0;
-    elmat(0,1)->SetSize(3*dof_u,dof_e); (*elmat(0,1))=0.0;
-    elmat(1,0)->SetSize(dof_e,3*dof_u); (*elmat(1,0))=0.0;
-    elmat(1,1)->SetSize(dof_e,dof_e); (*elmat(1,1))=0.0;
-
-
-    int dim =  Tr.GetSpaceDim();
-    if (dim != 3)
-    {
-        mfem::mfem_error("NLElPlastIntegrator::AssembleElementVector"
-                                " is not defined on manifold meshes");
-    }
-
-    Vector uu(elfun[0]->GetData()+0*dof_u, dof_u);
-    Vector vv(elfun[0]->GetData()+1*dof_u, dof_u);
-    Vector ww(elfun[0]->GetData()+2*dof_u, dof_u);
-
-    Vector ep(elfun[1]->GetData(), dof_e);
-
-    // temp storages for vectors and matrices
-    Vector su(dof_u); //shape functions for displacements
-    DenseMatrix du(dof_u,dim); //gradients of the shape functions
-    Vector dux; dux.SetDataAndSize(du.GetData()+0*dof_u,dof_u);
-    Vector duy; duy.SetDataAndSize(du.GetData()+1*dof_u,dof_u);
-    Vector duz; duz.SetDataAndSize(du.GetData()+2*dof_u,dof_u);
-    Vector se(dof_e); //shape functions for plastic strains
-    DenseMatrix de(dof_e,dim); //gradients of the shape functions
-    Vector grade(dim);
-
-    DenseMatrix vpss; //plastic strain at the integration points
-    DenseMatrix vkap; //accumulated plastic strain at the integration points
-    eep->GetElementValues(Tr.ElementNo,vpss);
-    kappa->GetElementValues(Tr.ElementNo,vkap);
-
-
-    DenseMatrix estrain(3,3);
-    Vector vestrain; vestrain.SetDataAndSize(estrain.GetData(),9);
-    DenseMatrix pstrain(3,3);
-    Vector vpstrain; vpstrain.SetDataAndSize(pstrain.GetData(),9);
-    DenseMatrix gradu(3,3); gradu=0.0;
-    Vector tv;
-
-    const IntegrationRule& ir=eep->GetElementIntRule(Tr.ElementNo);
-
-    double H=0.001;
-    double beta=0.0;
-    Vector vip(2); vip=0.0;
-
-    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> seval(&mat,&yf);
-    Vector stress(9);stress=0.0;
-
-    DenseMatrix Cep(9,9);
-    DenseMatrix tm(3*dof_u,3*dof_u);
-
-    DenseMatrix B(3*dof_u,9); B=0.0;
-    DenseMatrix Bm(3*dof_u,9);
-
-
-    double w;
-    for (int i = 0; i < ir.GetNPoints(); i++)
-    {
-        const IntegrationPoint &ip = ir.IntPoint(i);
-        Tr.SetIntPoint(&ip);
-        w = Tr.Weight();
-        w = ip.weight * w;
-
-        //calculate the total strain
-        el[0]->CalcPhysShape(Tr,su);
-        el[0]->CalcPhysDShape(Tr,du);
-        tv.SetDataAndSize(gradu.GetData()+0*3,dim);
-        du.MultTranspose(uu,tv);
-        tv.SetDataAndSize(gradu.GetData()+1*3,dim);
-        du.MultTranspose(vv,tv);
-        tv.SetDataAndSize(gradu.GetData()+2*3,dim);
-        du.MultTranspose(ww,tv);
-        //compute the strain tensor
-        for(int ii=0;ii<dim;ii++){
-            for(int jj=ii+1;jj<dim;jj++){
-                estrain(ii,jj)=0.5*(gradu(ii,jj)+gradu(jj,ii));
-                estrain(jj,ii)=estrain(ii,jj);
-            }
-            estrain(ii,ii)=gradu(ii,ii);
-        }
-
-        //estrain.PrintMatlab(std::cout);
-
-        //set current plastic strain Voight indexing
-        {
-            pstrain(0,0)=vpss(0,i);
-            pstrain(1,1)=vpss(1,i);
-            pstrain(2,2)=vpss(2,i);
-            pstrain(1,2)=pstrain(2,1)=vpss(3,i);
-            pstrain(0,2)=pstrain(2,0)=vpss(4,i);
-            pstrain(0,1)=pstrain(1,0)=vpss(5,i);
-        }
-
-        //calculate the filtered plastic strain
-        el[1]->CalcPhysShape(Tr,se);
-        el[1]->CalcPhysDShape(Tr,de);
-        double epf=se*ep;//filtered accumulated plastic strain
-        double rry=ll->Eval(Tr,ip); //filter radius
-
-        //evaluate the material behaviour at the integration point
-        {
-            double Em=E->Eval(Tr,ip);
-            double nm=nu->Eval(Tr,ip);
-            double ssy=ss_y->Eval(Tr,ip);
-            mat.SetE(Em);
-            mat.SetPoisson(nm);
-            yf.Set(ssy,H,beta);
-
-            vip[0]=vkap(0,i); //set the accumulated plastic starin
-            vip[1]=epf;   //set the filtered plastic strain
-            seval.SetStrain(vestrain);
-            seval.SetPlasticStrain(vpstrain);
-            seval.SetInternalParameters(vip);
-            //seval.Solve(stress,vpstrain,vip,Tr,ip);
-            seval.EvalTangent(Cep,Tr,ip);
-
-            //stress.Print(std::cout,3);
-            //Cep.PrintMatlab(std::cut);
-            //std::cout<<std::endl;
-
-        }
-
-        {
-            //form B
-            for(int i=0;i<dof_u;i++)
-            {
-                B(i+0*dof_u,0)=dux[i];
-                B(i+0*dof_u,1)=0.5*duy[i]; B(i+1*dof_u,1)=0.5*dux[i];
-                B(i+0*dof_u,2)=0.5*duz[i]; B(i+2*dof_u,2)=0.5*dux[i];
-                B(i+0*dof_u,3)=0.5*duy[i]; B(i+1*dof_u,3)=0.5*dux[i];
-                B(i+1*dof_u,4)=duy[i];
-                B(i+1*dof_u,5)=0.5*duz[i]; B(i+2*dof_u,5)=0.5*duy[i];
-                B(i+0*dof_u,6)=0.5*duz[i]; B(i+2*dof_u,6)=0.5*dux[i];
-                B(i+1*dof_u,7)=0.5*duz[i]; B(i+2*dof_u,7)=0.5*duy[i];
-                B(i+2*dof_u,8)=duz[i];
-            }
-
-            MultABt(B,Cep,Bm);
-            MultABt(B,Bm,tm);
-            elmat(0,0)->Add(w,tm);
-        }
-
-
-    }
-
-
-    //check the linear elasticity
-    /*
-    {
-        const IntegrationPoint &ip = ir.IntPoint(0);
-        double Em=E->Eval(Tr,ip);
-        double num=nu->Eval(Tr,ip);
-        double la=Em*num/((1.0+num)*(1.0-2.0*num));
-        double mu=Em/(2.0*(1+num));
-        ConstantCoefficient lc(la);
-        ConstantCoefficient mc(mu);
-        ElasticityIntegrator eli(lc,mc);
-        DenseMatrix K(dim*dof_u);
-        eli.AssembleElementMatrix(*el[0],Tr,K);
-        //std::fstream osi("k_mat.dat",std::ios::out);
-        //K.PrintMatlab(osi);
-
-        //K.Add(-1.0,*(elmat(0,0)));
-        //std::cout<<"|K|"<<K.FNorm()<<std::endl;
-        (*elmat(0,0))=K;
-
-    }
-    */
-
-
-    {
-        DenseMatrix K(dof_e,dof_e); K=0.0;
-        PowerCoefficient pc(*ll,2.0);
-        ConstantCoefficient one(1.0);
-        MassIntegrator mi(one);
-        mi.AssembleElementMatrix(*el[1],Tr,*elmat(1,1));
-        DiffusionIntegrator di(pc);
-        di.AssembleElementMatrix(*el[1],Tr,K);
-        elmat(1,1)->Add(1.0,K);
-    }
-
-}
-
-
-void NLElPlastIntegratorS::AssembleElementVector(const FiniteElement &el,
-                                                 ElementTransformation &Tr,
-                                                 const Vector &elfun, Vector &elvect)
-{
-    //the integrator works only for 3 dimensional problems
-    int dof_u = el.GetDof();
-    int dim =  Tr.GetSpaceDim();
-    elvect.SetSize(dim*dof_u); elvect=0.0;
-
-    Vector uu(elfun.GetData()+0*dof_u, dof_u);
-    Vector vv(elfun.GetData()+1*dof_u, dof_u);
-    Vector ww(elfun.GetData()+2*dof_u, dof_u);
-
-    Vector ru(elvect.GetData()+0*dof_u, dof_u); ru=0.0;
-    Vector rv(elvect.GetData()+1*dof_u, dof_u); rv=0.0;
-    Vector rw(elvect.GetData()+2*dof_u, dof_u); rw=0.0;
-
-    // temp storages for vectors and matrices
-    Vector su(dof_u); //shape functions for displacements
-    DenseMatrix du(dof_u,dim); //gradients of the shape functions
-    Vector dux; dux.SetDataAndSize(du.GetData()+0*dof_u,dof_u);
-    Vector duy; duy.SetDataAndSize(du.GetData()+1*dof_u,dof_u);
-    Vector duz; duz.SetDataAndSize(du.GetData()+2*dof_u,dof_u);
-
-    DenseMatrix vpss; //plastic strain at the integration points
-    DenseMatrix vkap; //accumulated plastic strain at the integration points
-    eep->GetElementValues(Tr.ElementNo,vpss);
-    kappa->GetElementValues(Tr.ElementNo,vkap);
-
-    DenseMatrix estrain(3,3);
-    Vector vestrain; vestrain.SetDataAndSize(estrain.GetData(),9);
-    DenseMatrix pstrain(3,3);
-    Vector vpstrain; vpstrain.SetDataAndSize(pstrain.GetData(),9);
-    DenseMatrix gradu(3,3); gradu=0.0;
-    Vector tv;
-
-    const IntegrationRule& ir=eep->GetElementIntRule(Tr.ElementNo);
-
-    double H=0.001;
-    double beta=0.0;
-    Vector vip(2); vip=0.0;
-
-    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> seval(&mat,&yf);
-    Vector stress(9);stress=0.0;
-
-    double w;
-    for (int i = 0; i < ir.GetNPoints(); i++)
-    {
-        const IntegrationPoint &ip = ir.IntPoint(i);
-        Tr.SetIntPoint(&ip);
-        w = Tr.Weight();
-        w = ip.weight * w;
-
-
-        //calculate the total strain
-        el.CalcPhysShape(Tr,su);
-        el.CalcPhysDShape(Tr,du);
-        tv.SetDataAndSize(gradu.GetData()+0*3,dim);
-        du.MultTranspose(uu,tv);
-        tv.SetDataAndSize(gradu.GetData()+1*3,dim);
-        du.MultTranspose(vv,tv);
-        tv.SetDataAndSize(gradu.GetData()+2*3,dim);
-        du.MultTranspose(ww,tv);
-        //compute the strain tensor
-        for(int ii=0;ii<dim;ii++){
-            for(int jj=ii+1;jj<dim;jj++){
-                estrain(ii,jj)=0.5*(gradu(ii,jj)+gradu(jj,ii));
-                estrain(jj,ii)=estrain(ii,jj);
-            }
-            estrain(ii,ii)=gradu(ii,ii);
-        }
-
-        //set current plastic strain Voight indexing
-        {
-            pstrain(0,0)=vpss(0,i);
-            pstrain(1,1)=vpss(1,i);
-            pstrain(2,2)=vpss(2,i);
-            pstrain(1,2)=pstrain(2,1)=vpss(3,i);
-            pstrain(0,2)=pstrain(2,0)=vpss(4,i);
-            pstrain(0,1)=pstrain(1,0)=vpss(5,i);
-        }
-
-        double epf=eef->GetValue(Tr,ip);
-
-        //evaluate the material behaviour at the integration point
-        {
-            double Em=E->Eval(Tr,ip);
-            double nm=nu->Eval(Tr,ip);
-            double ssy=ss_y->Eval(Tr,ip);
-            mat.SetE(Em);
-            mat.SetPoisson(nm);
-            yf.Set(ssy,H,beta);
-
-            vip[0]=vkap(0,i); //set the accumulated plastic starin
-            vip[1]=epf;   //set the filtered plastic strain
-            seval.SetStrain(vestrain);
-            seval.SetPlasticStrain(vpstrain);
-            seval.SetInternalParameters(vip);
-            seval.Solve(stress,vpstrain,vip,Tr,ip);
-            //stress.Print(std::cout,3);
-            //std::cout<<std::endl;
-        }
-
-        //assemble the RHS for displacements
-        ru.Add(stress[0]*w,dux);
-        ru.Add(0.5*stress[1]*w,duy); rv.Add(0.5*stress[1]*w,dux);
-        ru.Add(0.5*stress[2]*w,duz); rw.Add(0.5*stress[2]*w,dux);
-        ru.Add(0.5*stress[3]*w,duy); rv.Add(0.5*stress[3]*w,dux);
-        rv.Add(stress[4]*w,duy);
-        rv.Add(0.5*stress[5]*w,duz); rw.Add(0.5*stress[5]*w,duy);
-        ru.Add(0.5*stress[6]*w,duz); rw.Add(0.5*stress[6]*w,dux);
-        rv.Add(0.5*stress[7]*w,duz); rw.Add(0.5*stress[7]*w,duy);
-        rw.Add(stress[8]*w,duz);
-    }
-}
-
-void NLElPlastIntegratorS::AssembleElementGrad (const FiniteElement &el,
-                                                ElementTransformation &Tr,
-                                                const Vector &elfun, DenseMatrix &elmat)
-{
-    //the integrator works only for 3 dimensional problems
-    int dof_u = el.GetDof();
-    int dim =  Tr.GetSpaceDim();
-
-    elmat.SetSize(dof_u,dof_u); elmat=0.0;
-
-    if (dim != 3)
-    {
-        mfem::mfem_error("NLElPlastIntegrator::AssembleElementVector"
-                                " is not defined on manifold meshes");
-    }
-
-    Vector uu(elfun.GetData()+0*dof_u, dof_u);
-    Vector vv(elfun.GetData()+1*dof_u, dof_u);
-    Vector ww(elfun.GetData()+2*dof_u, dof_u);
-
-    // temp storages for vectors and matrices
-    Vector su(dof_u); //shape functions for displacements
-    DenseMatrix du(dof_u,dim); //gradients of the shape functions
-    Vector dux; dux.SetDataAndSize(du.GetData()+0*dof_u,dof_u);
-    Vector duy; duy.SetDataAndSize(du.GetData()+1*dof_u,dof_u);
-    Vector duz; duz.SetDataAndSize(du.GetData()+2*dof_u,dof_u);
-
-    DenseMatrix vpss; //plastic strain at the integration points
-    DenseMatrix vkap; //accumulated plastic strain at the integration points
-    eep->GetElementValues(Tr.ElementNo,vpss);
-    kappa->GetElementValues(Tr.ElementNo,vkap);
-
-
-    DenseMatrix estrain(3,3);
-    Vector vestrain; vestrain.SetDataAndSize(estrain.GetData(),9);
-    DenseMatrix pstrain(3,3);
-    Vector vpstrain; vpstrain.SetDataAndSize(pstrain.GetData(),9);
-    DenseMatrix gradu(3,3); gradu=0.0;
-    Vector tv;
-
-
-
-
-    const IntegrationRule& ir=eep->GetElementIntRule(Tr.ElementNo);
-
-
-    double H=0.001;
-    double beta=0.0;
-    Vector vip(2); vip=0.0;
-
-
-    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> seval(&mat,&yf);
-    Vector stress(9);stress=0.0;
-
-    DenseMatrix Cep(9,9);
-    DenseMatrix tm(3*dof_u,3*dof_u);
-
-    DenseMatrix B(3*dof_u,9); B=0.0;
-    DenseMatrix Bm(3*dof_u,9);
-
-    double w;
-    for (int i = 0; i < ir.GetNPoints(); i++)
-    {
-        const IntegrationPoint &ip = ir.IntPoint(i);
-        Tr.SetIntPoint(&ip);
-        w = Tr.Weight();
-        w = ip.weight * w;
-
-        //calculate the total strain
-        el.CalcPhysShape(Tr,su);
-        el.CalcPhysDShape(Tr,du);
-        tv.SetDataAndSize(gradu.GetData()+0*3,dim);
-        du.MultTranspose(uu,tv);
-        tv.SetDataAndSize(gradu.GetData()+1*3,dim);
-        du.MultTranspose(vv,tv);
-        tv.SetDataAndSize(gradu.GetData()+2*3,dim);
-        du.MultTranspose(ww,tv);
-        //compute the strain tensor
-        for(int ii=0;ii<dim;ii++){
-            for(int jj=ii+1;jj<dim;jj++){
-                estrain(ii,jj)=0.5*(gradu(ii,jj)+gradu(jj,ii));
-                estrain(jj,ii)=estrain(ii,jj);
-            }
-            estrain(ii,ii)=gradu(ii,ii);
-        }
-
-        //set current plastic strain Voight indexing
-        {
-            pstrain(0,0)=vpss(0,i);
-            pstrain(1,1)=vpss(1,i);
-            pstrain(2,2)=vpss(2,i);
-            pstrain(1,2)=pstrain(2,1)=vpss(3,i);
-            pstrain(0,2)=pstrain(2,0)=vpss(4,i);
-            pstrain(0,1)=pstrain(1,0)=vpss(5,i);
-        }
-
-        double epf=eef->GetValue(Tr,ip);
-        //evaluate the material behaviour at the integration point
-        {
-            double Em=E->Eval(Tr,ip);
-            double nm=nu->Eval(Tr,ip);
-            double ssy=ss_y->Eval(Tr,ip);
-            mat.SetE(Em);
-            mat.SetPoisson(nm);
-            yf.Set(ssy,H,beta);
-
-            vip[0]=vkap(0,i); //set the accumulated plastic starin
-            vip[1]=epf;   //set the filtered plastic strain
-            seval.SetStrain(vestrain);
-            seval.SetPlasticStrain(vpstrain);
-            seval.SetInternalParameters(vip);
-            //seval.Solve(stress,vpstrain,vip,Tr,ip);
-            seval.EvalTangent(Cep,Tr,ip);
-
-            //stress.Print(std::cout,3);
-            //Cep.PrintMatlab(std::cut);
-            //std::cout<<std::endl;
-
-        }
-
-        {
-            //form B
-            for(int i=0;i<dof_u;i++)
-            {
-                B(i+0*dof_u,0)=dux[i];
-                B(i+0*dof_u,1)=0.5*duy[i]; B(i+1*dof_u,1)=0.5*dux[i];
-                B(i+0*dof_u,2)=0.5*duz[i]; B(i+2*dof_u,2)=0.5*dux[i];
-                B(i+0*dof_u,3)=0.5*duy[i]; B(i+1*dof_u,3)=0.5*dux[i];
-                B(i+1*dof_u,4)=duy[i];
-                B(i+1*dof_u,5)=0.5*duz[i]; B(i+2*dof_u,5)=0.5*duy[i];
-                B(i+0*dof_u,6)=0.5*duz[i]; B(i+2*dof_u,6)=0.5*dux[i];
-                B(i+1*dof_u,7)=0.5*duz[i]; B(i+2*dof_u,7)=0.5*duy[i];
-                B(i+2*dof_u,8)=duz[i];
-            }
-
-            MultABt(B,Cep,Bm);
-            MultABt(B,Bm,tm);
-            elmat.Add(w,tm);
-        }
-
-    }
-
-}
-
-
-}

miniapps/elplast/elplast.hpp

@@ -1,310 +0,0 @@
-#ifndef ELPLAST_HPP
-#define ELPLAST_HPP
-
-#include "mfem.hpp"
-#include "coefficients.hpp"
-
-namespace mfem{
-
-
-
-
-class NLElPlastIntegrator:public BlockNonlinearFormIntegrator
-{
-public:
-    NLElPlastIntegrator(double E_=1, double nu_=0.2, double ss_y_=1.0, double ll_=0.1)
-    {
-        lE=new ConstantCoefficient(E_);
-        lnu=new ConstantCoefficient(nu_);
-        lss_y=new ConstantCoefficient(ss_y_);
-        lll=new ConstantCoefficient(ll_);
-
-        E=lE;
-        nu=lnu;
-        ss_y=lss_y;
-        ll=lll;
-
-        eep=nullptr;
-        kappa=nullptr;
-
-        force=nullptr;
-        flag_update=false;
-    }
-
-    virtual ~NLElPlastIntegrator()
-    {
-        delete lE;
-        delete lnu;
-        delete lss_y;
-        delete lll;
-    }
-
-    void SetPlasticStrains(mfem::QuadratureFunction& eep_,
-                           mfem::QuadratureFunction& kappa_)
-    {
-        eep=&eep_;
-        kappa=&kappa_;
-    }
-
-    void SetForce(VectorCoefficient& vc)
-    {
-        force=&vc;
-    }
-
-    void SetUpdateFlag(bool fl)
-    {
-        flag_update=fl;
-    }
-
-    virtual
-    double GetElementEnergy(const Array<const FiniteElement *> &el,
-                            ElementTransformation &Tr,
-                            const Array<const Vector *> &elfun)
-    {
-        return 0.0;
-    }
-
-
-    virtual
-    void AssembleElementVector(const Array<const FiniteElement *> &el,
-                               ElementTransformation &Tr,
-                               const Array<const Vector *> &elfun,
-                               const Array<Vector *> &elvec);
-
-    virtual
-    void AssembleElementGrad(const Array<const FiniteElement *> &el,
-                             ElementTransformation &Tr,
-                             const Array<const Vector *> &elfun,
-                             const Array2D<DenseMatrix *> &elmats);
-
-private:
-
-    mfem::Coefficient* lE;
-    mfem::Coefficient* lnu;
-    mfem::Coefficient* lss_y;
-    mfem::Coefficient* lll;
-
-
-    mfem::Coefficient* E;
-    mfem::Coefficient* nu;
-    mfem::Coefficient* ss_y;
-    mfem::Coefficient* ll;
-
-    mfem::VectorCoefficient* force;
-
-    mfem::QuadratureFunction* eep;
-    mfem::QuadratureFunction* kappa;
-
-    mfem::IsoElastMat mat;
-    mfem::J2YieldFunction yf;
-
-    bool flag_update;
-
-};
-
-class NLElPlastIntegratorS:public NonlinearFormIntegrator
-{
-public:
-
-    NLElPlastIntegratorS(double E_=1, double nu_=0.2, double ss_y_=1.0, double ll_=0.1)
-    {
-        lE=new ConstantCoefficient(E_);
-        lnu=new ConstantCoefficient(nu_);
-        lss_y=new ConstantCoefficient(ss_y_);
-        lll=new ConstantCoefficient(ll_);
-
-        E=lE;
-        nu=lnu;
-        ss_y=lss_y;
-        ll=lll;
-
-        eep=nullptr;
-        kappa=nullptr;
-    }
-
-    virtual ~NLElPlastIntegratorS()
-    {
-        delete lE;
-        delete lnu;
-        delete lss_y;
-        delete lll;
-    }
-
-    void SetPlasticStrains(mfem::QuadratureFunction& eep_,
-                           mfem::QuadratureFunction& kappa_)
-    {
-        eep=&eep_;
-        kappa=&kappa_;
-    }
-
-    void SetFilteredPlasticStrain(GridFunction& epf)
-    {
-        eef=&epf;
-    }
-
-    virtual
-    double GetElementEnrgy(const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
-    {
-        return 0.0;
-    }
-
-    virtual
-    void AssembleElementVector (const FiniteElement &el,
-                                ElementTransformation &Tr,
-                                const Vector &elfun, Vector &elvect);
-
-    virtual
-    void AssembleElementGrad (const FiniteElement &el,
-                              ElementTransformation &Tr,
-                              const Vector &elfun, DenseMatrix &elmat);
-
-
-private:
-    mfem::Coefficient* lE;
-    mfem::Coefficient* lnu;
-    mfem::Coefficient* lss_y;
-    mfem::Coefficient* lll;
-
-
-    mfem::Coefficient* E;
-    mfem::Coefficient* nu;
-    mfem::Coefficient* ss_y;
-    mfem::Coefficient* ll;
-
-    mfem::QuadratureFunction* eep;
-    mfem::QuadratureFunction* kappa;
-    mfem::GridFunction* eef;
-
-    mfem::IsoElastMat mat;
-    mfem::J2YieldFunction yf;
-};
-
-
-
-class ElPlastSolver
-{
-public:
-    ElPlastSolver(mfem::ParMesh* mesh_,int vorder=1, int forder=2);
-
-    ~ElPlastSolver();
-
-    /// Set the Newton Solver
-    void SetNewtonSolver(double rtol=1e-7, double atol=1e-12,int miter=1000, int prt_level=1);
-
-    /// Set the Linear Solver
-    void SetLinearSolver(double rtol=1e-8, double atol=1e-12, int miter=1000);
-
-    /// Solves the forward problem.
-    void FSolve();
-
-    /// Solves the adjoint with the provided rhs.
-    void ASolve(mfem::Vector& rhs);
-
-    /// Adds displacement BC in direction 0(x),1(y),2(z), or 4(all).
-    void AddDispBC(int id, int dir, double val);
-
-    /// Adds displacement BC in direction 0(x),1(y),2(z), or 4(all).
-    void AddDispBC(int id, int dir, mfem::Coefficient& val);
-
-    /// Adds displacement BC specified by the vector coefficient val.
-    void AddDispBC(int id, mfem::VectorCoefficient& val);
-
-    /// Adds vol force
-    void AddVolForce(int id, double fx, double fy, double fz);
-
-    /// Adds vol force
-    void AddVolForce(int id, mfem::VectorCoefficient& ff);
-
-    /// Returns the displacements.
-    mfem::ParGridFunction& GetDisplacements()
-    {
-        fdisp.SetFromTrueDofs(sol.GetBlock(0));
-        return fdisp;
-    }
-
-    /// Returns the adjoint displacements.
-    mfem::ParGridFunction& GetADisplacements()
-    {
-        adisp.SetFromTrueDofs(adj.GetBlock(0));
-        return adisp;
-    }
-
-    /// Returns the solution vector.
-    mfem::Vector& GetSol(){return sol;}
-
-    /// Returns the adjoint solution vector.
-    mfem::Vector& GetAdj(){return adj;}
-
-    void GetSol(ParGridFunction& sgf){
-        sgf.SetSpace(vfes); sgf.SetFromTrueDofs(sol.GetBlock(0));}
-
-    void GetAdj(ParGridFunction& agf){
-        agf.SetSpace(vfes); agf.SetFromTrueDofs(adj.GetBlock(0));}
-private:
-
-    double current_time;
-
-    mfem::ParMesh* pmesh;
-
-    //solution vector
-    mfem::BlockVector sol;
-    //adjoint vector
-    mfem::BlockVector adj;
-    //RHS
-    mfem::BlockVector rhs;
-
-    // localy defined volumetric forces
-    std::map<int, mfem::VectorConstantCoefficient*> lvforce;
-    // globaly defined volumetric forces
-    std::map<int, mfem::VectorCoefficient*> volforce;
-
-    // boundary conditions for x,y, and z directions
-    std::map<int, mfem::ConstantCoefficient> bcx;
-    std::map<int, mfem::ConstantCoefficient> bcy;
-    std::map<int, mfem::ConstantCoefficient> bcz;
-
-    // holds BC in coefficient form
-    std::map<int, mfem::Coefficient*> bccx;
-    std::map<int, mfem::Coefficient*> bccy;
-    std::map<int, mfem::Coefficient*> bccz;
-    std::map<int, mfem::VectorCoefficient*> bcca;
-
-    // holds the displacement contrained DOFs
-    mfem::Array<int> ess_tdofv;
-
-    //forward solution
-    mfem::ParGridFunction fdisp;
-    //adjoint solution
-    mfem::ParGridFunction adisp;
-
-    //total strains
-    mfem::QuadratureFunction eee;
-    //plastic strains
-    mfem::QuadratureFunction eep;
-    //accumulated plastic strain
-    mfem::QuadratureFunction kappa;
-
-    //Newton solver parameters
-    double abs_tol;
-    double rel_tol;
-    int print_level;
-    int max_iter;
-
-    //Linear solver parameters
-    double linear_rtol;
-    double linear_atol;
-    int linear_iter;
-
-    mfem::ParBlockNonlinearForm *nf;
-    mfem::ParFiniteElementSpace* vfes; //displacements fes
-    mfem::ParFiniteElementSpace* ffes; //filter fes
-
-    mfem::FiniteElementCollection* vfec;
-    mfem::FiniteElementCollection* ffec;
-
-    mfem::QuadratureSpace* qfes;
-};
-
-}
-
-#endif

miniapps/elplast/holep.cpp

@@ -1,401 +0,0 @@
-//                       MFEM Example 21 - Parallel Version
-//
-// Compile with: make ex21p
-//
-// Sample runs:  mpirun -np 4 ex21p
-//               mpirun -np 4 ex21p -o 3
-//               mpirun -np 4 ex21p -m ../data/beam-quad.mesh
-//               mpirun -np 4 ex21p -m ../data/beam-quad.mesh -o 3
-//               mpirun -np 4 ex21p -m ../data/beam-tet.mesh
-//               mpirun -np 4 ex21p -m ../data/beam-tet.mesh -o 2
-//               mpirun -np 4 ex21p -m ../data/beam-hex.mesh
-//               mpirun -np 4 ex21p -m ../data/beam-hex.mesh -o 2
-//
-// Description:  This is a version of Example 2p with a simple adaptive mesh
-//               refinement loop. The problem being solved is again the linear
-//               elasticity describing a multi-material cantilever beam.
-//               The problem is solved on a sequence of meshes which
-//               are locally refined in a conforming (triangles, tetrahedrons)
-//               or non-conforming (quadrilaterals, hexahedra) manner according
-//               to a simple ZZ error estimator.
-//
-//               The example demonstrates MFEM's capability to work with both
-//               conforming and nonconforming refinements, in 2D and 3D, on
-//               linear and curved meshes. Interpolation of functions from
-//               coarse to fine meshes, as well as persistent GLVis
-//               visualization are also illustrated.
-//
-//               We recommend viewing Examples 2p and 6p before viewing this
-//               example.
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   // 0. Initialize MPI and HYPRE.
-   Mpi::Init(argc, argv);
-   int num_procs = Mpi::WorldSize();
-   int myid = Mpi::WorldRank();
-   Hypre::Init();
-
-   // 1. Parse command-line options.
-   const char *mesh_file = "hole-abaqus.mesh";
-   int serial_ref_levels = 0;
-   int order = 2;
-   //bool static_cond = false;
-   bool static_cond = true;
-   bool visualization = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&serial_ref_levels, "-rs", "--refine-serial",
-                  "Number of uniform serial refinements (before parallel"
-                  " partitioning)");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree).");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   // 2. Read the mesh from the given mesh file. We can handle triangular,
-   //    quadrilateral, tetrahedral, and hexahedral meshes with the same code.
-   Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-   MFEM_VERIFY(mesh.SpaceDimension() == dim, "invalid mesh");
-
-   std::cout<<"NE="<<mesh.GetNE()<<std::endl;
-   std::cout<<"NN="<<mesh.GetNV()<<std::endl;
-
-   double* v1=mesh.GetVertex(0);
-   std::cout<<v1[0]<<" "<<v1[1]<<" "<<v1[2]<<std::endl;
-
-
-   // 3. Refine the mesh before parallel partitioning. Since a NURBS mesh can
-   //    currently only be refined uniformly, we need to convert it to a
-   //    piecewise-polynomial curved mesh. First we refine the NURBS mesh a bit
-   //    more and then project the curvature to quadratic Nodes.
-   if (mesh.NURBSext && serial_ref_levels == 0)
-   {
-      serial_ref_levels = 2;
-   }
-   for (int i = 0; i < serial_ref_levels; i++)
-   {
-      mesh.UniformRefinement();
-   }
-   if (mesh.NURBSext)
-   {
-      mesh.SetCurvature(2);
-   }
-   mesh.EnsureNCMesh();
-
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   mesh.Clear();
-
-   // 4. Define a finite element space on the mesh. The polynomial order is
-   //    one (linear) by default, but this can be changed on the command line.
-   H1_FECollection fec(order, dim);
-   ParFiniteElementSpace fespace(&pmesh, &fec, dim);
-
-   // 5. As in Example 2, we set up the linear form b(.) which corresponds to
-   //    the right-hand side of the FEM linear system. In this case, b_i equals
-   //    the boundary integral of f*phi_i where f represents a "pull down"
-   //    force on the Neumann part of the boundary and phi_i are the basis
-   //    functions in the finite element fespace. The force is defined by the
-   //    VectorArrayCoefficient object f, which is a vector of Coefficient
-   //    objects. The fact that f is non-zero on boundary attribute 2 is
-   //    indicated by the use of piece-wise constants coefficient for its last
-   //    component. We don't assemble the discrete problem yet, this will be
-   //    done in the main loop.
-   VectorArrayCoefficient f(dim);
-   for (int i = 0; i < dim-1; i++)
-   {
-      f.Set(i, new ConstantCoefficient(0.0));
-   }
-
-      Vector pull_force(pmesh.attributes.Max());
-      pull_force = 0.0;
-   // pull_force(1) = -1.0e-2;
-      pull_force(1) = 0.757576E-2;
-      pull_force(2) = 0.757576E-2;
-
-      pull_force*=100000;
-
-      f.Set(0, new PWConstCoefficient(pull_force));
-
-   ParLinearForm b(&fespace);
-
-   b.AddDomainIntegrator(new VectorDomainLFIntegrator(f));
-
-
-   // 6. Set up the bilinear form a(.,.) on the finite element space
-   //    corresponding to the linear elasticity integrator with piece-wise
-   //    constants coefficient lambda and mu.
-
-   std::cout<<"pmesh atr="<<pmesh.attributes.Max()<<std::endl;
-
-   Vector lambda(pmesh.attributes.Max());
-   lambda = 1.0;
-   //lambda(0) = lambda(1)*50;
-   lambda(0) = 0.138651E+8;
-   lambda(1) = 0.576923E+7;
-   lambda(2) = 0.288462E+6;
-
-   PWConstCoefficient lambda_func(lambda);
-   Vector mu(pmesh.attributes.Max());
-   mu = 1.0;
-   //mu(0) = mu(1)*50;
-   mu(0) = 0.11811E+8;
-   mu(1) = 0.384615E+7;
-   mu(2) = 0.192308E+6;
-
-   PWConstCoefficient mu_func(mu);
-
-   ParBilinearForm a(&fespace);
-   BilinearFormIntegrator *integ =
-      new ElasticityIntegrator(lambda_func,mu_func);
-   a.AddDomainIntegrator(integ);
-   if (static_cond) { a.EnableStaticCondensation(); }
-
-   // 7. The solution vector x and the associated finite element grid function
-   //    will be maintained over the AMR iterations. We initialize it to zero.
-   Vector zero_vec(dim);
-   zero_vec = 0.0;
-   VectorConstantCoefficient zero_vec_coeff(zero_vec);
-   ParGridFunction x(&fespace);
-   x = 0.0;
-
-   // 8. Determine the list of true (i.e. conforming) essential boundary dofs.
-   //    In this example, the boundary conditions are defined by marking only
-   //    boundary attribute 1 from the mesh as essential and converting it to a
-   //    list of true dofs.  The conversion to true dofs will be done in the
-   //    main loop.
-   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
-   ess_bdr = 0;
-   ess_bdr[0] = 1;
-
-   // 9. GLVis visualization.
-   char vishost[] = "localhost";
-   int  visport   = 19916;
-   socketstream sol_sock;
-
-   // 10. Set up an error estimator. Here we use the Zienkiewicz-Zhu estimator
-   //     that uses the ComputeElementFlux method of the ElasticityIntegrator to
-   //     recover a smoothed flux (stress) that is subtracted from the element
-   //     flux to get an error indicator. We need to supply the space for the
-   //     smoothed flux: an (H1)^tdim (i.e., vector-valued) space is used here.
-   //     Here, tdim represents the number of components for a symmetric (dim x
-   //     dim) tensor.
-   const int tdim = dim*(dim+1)/2;
-   L2_FECollection flux_fec(order, dim);
-   ParFiniteElementSpace flux_fespace(&pmesh, &flux_fec, tdim);
-   ParFiniteElementSpace smooth_flux_fespace(&pmesh, &fec, tdim);
-   L2ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace,
-                                       smooth_flux_fespace);
-
-
-   // ParaView output.
-   ParaViewDataCollection dacol("ParaView", &pmesh);
-   dacol.SetLevelsOfDetail(order);
-   dacol.SetDataFormat(VTKFormat::ASCII);
-   dacol.RegisterField("disp", &x);
-   dacol.SetTime(1.0);
-   dacol.SetCycle(1);
-   dacol.Save();
-
-
-   // 11. A refiner selects and refines elements based on a refinement strategy.
-   //     The strategy here is to refine elements with errors larger than a
-   //     fraction of the maximum element error. Other strategies are possible.
-   //     The refiner will call the given error estimator.
-   ThresholdRefiner refiner(estimator);
-   refiner.SetTotalErrorFraction(0.7);
-
-   // 12. The main AMR loop. In each iteration we solve the problem on the
-   //     current mesh, visualize the solution, and refine the mesh.
-   const int max_dofs = 500000;
-   const int max_amr_itr = 20;
-   for (int it = 0; it <= max_amr_itr; it++)
-   {
-      HYPRE_BigInt global_dofs = fespace.GlobalTrueVSize();
-      if (myid == 0)
-      {
-         cout << "\nAMR iteration " << it << endl;
-         cout << "Number of unknowns: " << global_dofs << endl;
-      }
-
-      // 13. Assemble the stiffness matrix and the right-hand side.
-      a.Assemble();
-      b.Assemble();
-
-      // 14. Set Dirichlet boundary values in the GridFunction x.
-      //     Determine the list of Dirichlet true DOFs in the linear system.
-      Array<int> ess_tdof_list;
-      x.ProjectBdrCoefficient(zero_vec_coeff, ess_bdr);
-      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
-
-      // 15. Create the linear system: eliminate boundary conditions, constrain
-      //     hanging nodes and possibly apply other transformations. The system
-      //     will be solved for true (unconstrained) DOFs only.
-
-      HypreParMatrix A;
-      Vector B, X;
-      const int copy_interior = 1;
-      a.FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
-
-      // 16. Define and apply a parallel PCG solver for AX=B with the BoomerAMG
-      //     preconditioner from hypre.
-      HypreBoomerAMG amg;
-      amg.SetPrintLevel(0);
-      // amg.SetSystemsOptions(dim); // optional
-      CGSolver pcg(A.GetComm());
-      pcg.SetPreconditioner(amg);
-      pcg.SetOperator(A);
-      pcg.SetRelTol(1e-6);
-      pcg.SetMaxIter(5000);
-      pcg.SetPrintLevel(3); // print the first and the last iterations only
-      pcg.Mult(B, X);
-
-      // 17. After solving the linear system, reconstruct the solution as a
-      //     finite element GridFunction. Constrained nodes are interpolated
-      //     from true DOFs (it may therefore happen that x.Size() >= X.Size()).
-      a.RecoverFEMSolution(X, b, x);
-
-      //ParaView output
-      dacol.SetTime(double(it+1));
-      dacol.SetCycle(it+1);
-      dacol.Save();
-
-
-      // 18. Send solution by socket to the GLVis server.
-      if (visualization && it == 0)
-      {
-         sol_sock.open(vishost, visport);
-         sol_sock.precision(8);
-      }
-      if (visualization && sol_sock.good())
-      {
-         GridFunction nodes(&fespace), *nodes_p = &nodes;
-         pmesh.GetNodes(nodes);
-         nodes += x;
-         int own_nodes = 0;
-         pmesh.SwapNodes(nodes_p, own_nodes);
-         x.Neg(); // visualize the backward displacement
-         sol_sock << "parallel " << num_procs << ' ' << myid << '\n';
-         sol_sock << "solution\n" << pmesh << x << flush;
-         x.Neg();
-         pmesh.SwapNodes(nodes_p, own_nodes);
-         if (it == 0)
-         {
-            sol_sock << "keys '" << ((dim == 2) ? "Rjl" : "") << "m'" << endl;
-         }
-         sol_sock << "window_title 'AMR iteration: " << it << "'\n"
-                  << "pause" << endl;
-         if (myid == 0)
-         {
-            cout << "Visualization paused. "
-                 "Press <space> in the GLVis window to continue." << endl;
-         }
-      }
-
-      if (global_dofs > max_dofs)
-      {
-         if (myid == 0)
-         {
-            cout << "Reached the maximum number of dofs. Stop." << endl;
-         }
-         break;
-      }
-
-      // 19. Call the refiner to modify the mesh. The refiner calls the error
-      //     estimator to obtain element errors, then it selects elements to be
-      //     refined and finally it modifies the mesh. The Stop() method can be
-      //     used to determine if a stopping criterion was met.
-      refiner.Apply(pmesh);
-      if (refiner.Stop())
-      {
-         if (myid == 0)
-         {
-            cout << "Stopping criterion satisfied. Stop." << endl;
-         }
-         break;
-      }
-
-      // 20. Update the space to reflect the new state of the mesh. Also,
-      //     interpolate the solution x so that it lies in the new space but
-      //     represents the same function. This saves solver iterations later
-      //     since we'll have a good initial guess of x in the next step.
-      //     Internally, FiniteElementSpace::Update() calculates an
-      //     interpolation matrix which is then used by GridFunction::Update().
-      fespace.Update();
-      x.Update();
-
-      // 21. Load balance the mesh, and update the space and solution. Currently
-      //     available only for nonconforming meshes.
-      if (pmesh.Nonconforming())
-      {
-         pmesh.Rebalance();
-
-         // Update the space and the GridFunction. This time the update matrix
-         // redistributes the GridFunction among the processors.
-         fespace.Update();
-         x.Update();
-      }
-
-      // 21. Inform also the bilinear and linear forms that the space has
-      //     changed.
-      a.Update();
-      b.Update();
-   }
-
-   dacol.SetTime(double(max_amr_itr));
-   dacol.SetCycle(max_amr_itr);
-   dacol.Save();
-
-   {
-      ostringstream mref_name, mesh_name, sol_name;
-      mref_name << "ex21p_reference_mesh." << setfill('0') << setw(6) << myid;
-      mesh_name << "ex21p_deformed_mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "ex21p_displacement." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ref_out(mref_name.str().c_str());
-      mesh_ref_out.precision(16);
-      pmesh.Print(mesh_ref_out);
-
-      ofstream mesh_out(mesh_name.str().c_str());
-      mesh_out.precision(16);
-      GridFunction nodes(&fespace), *nodes_p = &nodes;
-      pmesh.GetNodes(nodes);
-      nodes += x;
-      int own_nodes = 0;
-      pmesh.SwapNodes(nodes_p, own_nodes);
-      pmesh.Print(mesh_out);
-      pmesh.SwapNodes(nodes_p, own_nodes);
-
-      ofstream x_out(sol_name.str().c_str());
-      x_out.precision(16);
-      x.Save(x_out);
-   }
-
-   return 0;
-}

miniapps/elplast/makefile

@@ -1,89 +0,0 @@
-# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-# Use the MFEM build directory
-MFEM_DIR ?= ../..
-MFEM_BUILD_DIR ?= ../..
-SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/elplast/,)
-CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
-# Use the MFEM install directory
-# MFEM_INSTALL_DIR = ../../mfem
-# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
-
-# Include defaults.mk to get XLINKER
-DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
-include $(DEFAULTS_MK)
-
-MFEM_LIB_FILE = mfem_is_not_built
--include $(CONFIG_MK)
-
-ADAPT_SRC = adaptive_el.cpp
-ADAPT_OBJ = $(ADAPT_SRC:.cpp=.o)
-
-PAR_MINIAPPS = adaptel
-
-ifeq ($(MFEM_USE_MPI),NO)
-   MINIAPPS = 
-else
-   MINIAPPS = $(PAR_MINIAPPS)
-endif
-
-COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
-
-# If MFEM_SHARED is set, add the ../common rpath
-COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
-
-.SUFFIXES:
-.SUFFIXES: .o .cpp .mk
-.PHONY: all lib-common clean clean-build clean-exec
-
-
-# Remove built-in rules
-%: %.cpp
-%.o: %.cpp
-
-%: %.o $(NAVIER_COMMON_OBJ)
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $^ -o $@ $(MFEM_LIBS)
-
-%.o: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_FLAGS) -c $< -o $@
-
-
-
-
-all: $(MINIAPPS)
-
-adaptel: $(ADAPT_OBJ) lib-common
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $(ADAPT_OBJ) $(COMMON_LIB) $(MFEM_LIBS)
-
-# Rule for building lib-common
-lib-common:
-	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
-
-MFEM_TESTS = MINIAPPS
-include $(MFEM_TEST_MK)
-
-
-# Generate an error message if the MFEM library is not built and exit
-$(MFEM_LIB_FILE):
-	$(error The MFEM library is not built)
-
-clean: clean-build clean-exec
-
-clean-build:
-	rm -f *.o *~ adaptel
-	rm -rf *.dSYM *.TVD.*breakpoints
-
-clean-exec:
-	@rm -f diffusion.mesh diffusion.gf
-	@rm -rf ParaViewDistance ParaViewDiffusion ParaViewExtrapolate ParaViewLSF

miniapps/elplast/test_elinteg.cpp

@@ -1,148 +0,0 @@
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-#include "coefficients.hpp"
-
-
-int main(int argc, char *argv[])
-{
-
-    mfem::ElementTransformation* Tr;
-    mfem::IntegrationPoint* ip;
-
-
-    mfem::IsoElastMat mat; mat.SetElastParam(1,0.0);
-    mfem::J2YieldFunction yf(1.0,0.10,1.0);
-
-    mfem::StressEval<mfem::IsoElastMat,mfem::J2YieldFunction> se(&mat,&yf);
-
-    mfem::Vector ss(9); ss=0.0;
-    mfem::Vector ep(9); ep=0.0;
-    mfem::Vector iv(2); iv=0.0;
-
-    mfem::Vector dee(9); dee=0.0;
-    mfem::Vector ee(9); ee=0.0;
-    mfem::Vector epn(9); epn=0.0;
-    mfem::Vector ivn(2); ivn=0.0;
-
-    int bb=8;
-
-    for(int i=0;i<20;i++){
-
-
-        std::cout<<i<<" " ;//std::endl;
-
-        dee(bb)=0.11;
-        ee.Add(1.0,dee);
-        se.SetStrain(ee);
-        se.SetPlasticStrain(epn);
-        se.SetInternalParameters(ivn);
-        se.Solve(ss,epn,ivn,*Tr,*ip);
-        //epn=ep;
-        //ivn=iv;
-
-        std::cout<<ee[bb]<<" "<<ss[bb]<<" "<<ss[1]<<" "<<ep[bb]<<" "<<yf.Eval(ss,iv)<<std::endl;
-
-        /*
-        std::cout<<"strain="<<std::endl;
-        ee.Print(std::cout,3);
-
-        std::cout<<"stress="<<std::endl;
-        std::cout.precision(6);
-        ss.Print(std::cout,3);
-
-        iv.Print(std::cout);
-        */
-
-    }
-
-    for(int i=0;i<40;i++){
-
-
-        std::cout<<i<<" " ;//std::endl;
-
-        dee(bb)=-0.11;
-        ee.Add(1.0,dee);
-        se.SetStrain(ee);
-        se.SetPlasticStrain(epn);
-        se.SetInternalParameters(ivn);
-        se.Solve(ss,ep,iv,*Tr,*ip);
-        epn=ep;
-        ivn=iv;
-
-        std::cout<<ee[bb]<<" "<<ss[bb]<<" "<<ss[1]<<" "<<ep[bb]<<" "<<yf.Eval(ss,iv)<<std::endl;
-
-        /*
-        std::cout<<"strain="<<std::endl;
-        ee.Print(std::cout,3);
-
-        std::cout<<"stress="<<std::endl;
-        std::cout.precision(6);
-        ss.Print(std::cout,3);
-
-        iv.Print(std::cout);
-        */
-
-    }
-
-    for(int i=0;i<25;i++){
-
-
-        std::cout<<i<<" " ;//std::endl;
-
-        dee(bb)=0.11;
-        ee.Add(1.0,dee);
-        se.SetStrain(ee);
-        se.SetPlasticStrain(epn);
-        se.SetInternalParameters(ivn);
-        se.Solve(ss,ep,iv,*Tr,*ip);
-        epn=ep;
-        ivn=iv;
-
-        std::cout<<ee[bb]<<" "<<ss[bb]<<" "<<ss[1]<<" "<<ep[bb]<<" "<<yf.Eval(ss,iv)<<std::endl;
-
-        /*
-        std::cout<<"strain="<<std::endl;
-        ee.Print(std::cout,3);
-
-        std::cout<<"stress="<<std::endl;
-        std::cout.precision(6);
-        ss.Print(std::cout,3);
-
-        iv.Print(std::cout);
-        */
-
-    }
-
-    for(int i=0;i<15;i++){
-
-
-        std::cout<<i<<" " ;//std::endl;
-
-        dee(bb)=-0.11;
-        ee.Add(1.0,dee);
-        se.SetStrain(ee);
-        se.SetPlasticStrain(epn);
-        se.SetInternalParameters(ivn);
-        se.Solve(ss,ep,iv,*Tr,*ip);
-        epn=ep;
-        ivn=iv;
-
-        std::cout<<ee[bb]<<" "<<ss[bb]<<" "<<ss[1]<<" "<<ep[bb]<<" "<<yf.Eval(ss,iv)<<std::endl;
-
-        /*
-        std::cout<<"strain="<<std::endl;
-        ee.Print(std::cout,3);
-
-        std::cout<<"stress="<<std::endl;
-        std::cout.precision(6);
-        ss.Print(std::cout,3);
-
-        iv.Print(std::cout);
-        */
-
-    }
-
-
-    return 0;
-}