Update bazel build modes

VectorDivergenceIntegrator on surface meshes

.bazelignore

@@ -0,0 +1,4 @@
+bazel-bazel
+bazel-bin
+bazel-out
+bazel-testlogs

.bazelrc

@@ -0,0 +1,44 @@
+# Bazel configuration file
+# https://docs.bazel.build/versions/main/guide.html#bazelrc
+# --compilation_mode [-c] (fastbuild, dbg or opt; default: "fastbuild")
+# --[no]subcommands [-s] (true, pretty_print or false; default: "false")
+# bazel clean --expunge
+
+# Enable Bzlmod for every Bazel command
+common --enable_bzlmod
+
+# Build options ###############################################################
+build --color=yes
+
+# Disk cache setup ############################################################
+build --disk_cache=/tmp/bazel_disk_cache
+
+# Alias definition ############################################################
+# Allow  to specify the build mode and precision with the following flags:
+# --mode=[serial|parallel]
+# --precision=[single|double]
+build --flag_alias=mode=//:mode
+build --flag_alias=precision=//:precision
+
+# Compilation options #########################################################
+build --copt=-g
+build --copt=-Wall
+# build --copt=-Wextra
+# build --copt=-Werror
+build --copt=-std=c++17
+build --cxxopt=-Wno-unused-parameter
+build --strip=never
+
+# build:macos --copt=-Werror
+# build:macos --cxxopt=-Wno-error=deprecated-declarations
+# build:macos --cxxopt=-Wno-error=deprecated-enum-enum-conversion
+# build:macos --cxxopt=-Wno-error=deprecated-anon-enum-enum-conversion
+# build:macos --cxxopt=-Wno-error=inconsistent-missing-override
+# build:macos --cxxopt=-Wno-error=overloaded-virtual
+
+# Prevent Bazel from detecting the system's C++ toolchain.
+# build --action_env=BAZEL_DO_NOT_DETECT_CPP_TOOLCHAIN=1
+# build --incompatible_strict_action_env=true
+
+# Enable the CC toolchain resolution based on platforms.
+# build --incompatible_enable_cc_toolchain_resolution

.gitignore

@@ -19,6 +19,9 @@ CMakeFiles/
 # Clangd server cache
 *.cache*
 
+# VSCode files
+.vscode/
+
 # Backup files
 *~
 
@@ -51,6 +54,11 @@ doc/warnings.log
 *.dSYM
 .DS_Store
 
+# Bazel specific
+bazel-*
+*.bazel_cache
+MODULE.bazel.lock
+
 # Example and miniapp binaries and outputs
 
 examples/ex[0-9]

BUILD

@@ -0,0 +1,468 @@
+### String Flags ##############################################################
+load("@bazel_skylib//rules:common_settings.bzl", "string_flag")
+
+### Load rules ################################################################
+load("@rules_cc//cc:defs.bzl", "cc_library")
+
+### Load config/bazel/config ##################################################
+load("//config/bazel:config.bzl", "mfem_serial_examples",
+                                  "mfem_parallel_examples", "mfem_use")
+
+### Load config/bazel/settings ################################################
+### https://bazel.build/docs/configurable-attributes
+load("//config/bazel:settings.bzl", "mode", "precision", "print_mode", "print_precision")
+
+# FLAG: Serial/Parallel MODE ##################################################
+string_flag(
+    name = "mode",
+    values = ["serial", "parallel"],
+    build_setting_default = "serial",
+)
+
+mode(name = "serial")
+
+mode(name = "parallel")
+
+config_setting(
+    name = "serial_mode",
+    flag_values = {":mode": "serial"},
+)
+
+config_setting(
+    name = "parallel_mode",
+    flag_values = {":mode": "parallel"},
+)
+
+print_mode(
+    name = "print_mode",
+    mode = select({
+        ":serial_mode": "serial",
+        ":parallel_mode": "parallel",
+    }),
+)
+
+# FLAG: Double/Single PRECISION ###############################################
+string_flag(
+    name = "precision",
+    values = ["single", "double"],
+    build_setting_default = "double",
+)
+
+precision(name = "double")
+
+precision(name = "single")
+
+config_setting(
+    name = "single_precision",
+    flag_values = {":precision": "single"},
+)
+
+config_setting(
+    name = "double_precision",
+    flag_values = {":precision": "double"},
+)
+
+print_precision(
+    name = "print_precision",
+    precision = select({
+        ":single_precision": "single",
+        ":double_precision": "double",
+    }),
+)
+
+# MFEM_USE_* definitions ######################################################
+mfem_use(
+    name = "mfem_not_mpi",
+    define = "MFEM_USE_MPI",
+    use = False,
+)
+
+mfem_use(
+    name = "mfem_use_mpi",
+    define = "MFEM_USE_MPI",
+    use = True,
+)
+
+mfem_use(
+    name = "mfem_not_metis",
+    define = "MFEM_USE_METIS",
+    use = False,
+)
+
+mfem_use(
+    name = "mfem_use_metis",
+    define = "MFEM_USE_METIS",
+    use = True,
+)
+
+mfem_use(
+    name = "mfem_not_metis_5",
+    define = "MFEM_USE_METIS_5",
+    use = False,
+)
+
+mfem_use(
+    name = "mfem_use_metis_5",
+    define = "MFEM_USE_METIS_5",
+    use = True,
+)
+
+mfem_use(
+    name = "no_mfem_hypre_version",
+    define = "MFEM_HYPRE_VERSION",
+    use = False,
+)
+
+mfem_use(
+    name = "mfem_hypre_version",
+    define = "MFEM_HYPRE_VERSION",
+    use = True,
+)
+
+mfem_use(
+    name = "mfem_not_double",
+    define = "MFEM_USE_DOUBLE",
+    use = False,
+)
+
+mfem_use(
+    name = "mfem_use_double",
+    define = "MFEM_USE_DOUBLE",
+    use = True,
+)
+
+mfem_use(
+    name = "mfem_not_single",
+    define = "MFEM_USE_SINGLE",
+    use = False,
+)
+
+mfem_use(
+    name = "mfem_use_single",
+    define = "MFEM_USE_SINGLE",
+    use = True,
+)
+
+### https://bazel.build/reference/be/general#genrule
+genrule(
+    name = "genrule_config_bazel",
+    srcs = ["BUILD"],
+    outs = ["config/bazel.hpp"],
+    cmd = """cat <<EOF > $@
+// Copyright (c) 2010-2025, 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.
+
+#ifndef MFEM_CONFIG_HEADER
+#define MFEM_CONFIG_HEADER
+
+// MFEM version: integer of the form: (major*100 + minor)*100 + patch.
+#define MFEM_VERSION 40701
+
+// MFEM version string of the form "3.3" or "3.3.1".
+#define MFEM_VERSION_STRING "4.7.1"
+
+// MFEM version type, see the MFEM_VERSION_TYPE_* constants below.
+#define MFEM_VERSION_TYPE ((MFEM_VERSION) % 2)
+
+// MFEM version type constants.
+#define MFEM_VERSION_TYPE_RELEASE 0
+#define MFEM_VERSION_TYPE_DEVELOPMENT 1
+
+// Separate MFEM version numbers for major, minor, and patch.
+#define MFEM_VERSION_MAJOR ((MFEM_VERSION) / 10000)
+#define MFEM_VERSION_MINOR (((MFEM_VERSION) / 100) % 100)
+#define MFEM_VERSION_PATCH ((MFEM_VERSION) % 100)
+
+// The absolute path of the MFEM source prefix.
+#define MFEM_SOURCE_DIR "$$(realpath $$(realpath BUILD)/..)"
+
+// The absolute path of the MFEM installation prefix.
+#define MFEM_INSTALL_DIR "$$(realpath $(BINDIR))"
+
+// Description of the git commit used to build MFEM.
+#define MFEM_GIT_STRING "heads/bazel-git-..."
+
+// Build the parallel MFEM library.
+// Requires an MPI compiler, and the libraries HYPRE and METIS.
+$(MFEM_USE_MPI)
+
+// Enable MFEM features that use the METIS library (parallel MFEM).
+$(MFEM_USE_METIS)
+
+// Enable this option if linking with METIS version 5 (parallel MFEM).
+$(MFEM_USE_METIS_5)
+
+// Version of HYPRE used for building MFEM.
+// macOS: 23200, ubuntu: 21821
+$(MFEM_HYPRE_VERSION) 23200
+
+// Use single/double-precision floating point type
+$(MFEM_USE_DOUBLE)
+$(MFEM_USE_SINGLE)
+
+// Internal MFEM option: enable group/batch allocation for some small objects.
+#define MFEM_USE_MEMALLOC
+
+// Which library functions to use in class StopWatch for measuring time.
+// For a list of the available options, see INSTALL.
+// If not defined, an option is selected automatically.
+// 0/1/2/3/4/5/6/NO
+#define MFEM_TIMER_TYPE 0
+
+#endif // MFEM_CONFIG_HEADER
+EOF""",
+    local = False,
+    message = "Generating config bazel.hpp file",
+    toolchains = select({
+        ":serial_mode": [
+            ":mfem_not_mpi",
+            ":mfem_not_metis",
+            ":mfem_not_metis_5",
+            ":no_mfem_hypre_version",
+        ],
+        ":parallel_mode": [
+            ":mfem_use_mpi",
+            ":mfem_use_metis",
+            ":mfem_use_metis_5",
+            ":mfem_hypre_version",
+        ],
+        "//conditions:default": [
+            ":mfem_not_mpi",
+            ":mfem_not_metis",
+            ":mfem_not_metis_5",
+            ":mfem_hypre_version",
+            ":mfem_use_double",
+            ":mfem_not_single",
+        ],
+    }) + select({
+        ":double_precision": [
+            ":mfem_use_double",
+            ":mfem_not_single",
+        ],
+        ":single_precision": [
+            ":mfem_not_double",
+            ":mfem_use_single",
+        ],
+    }),
+)
+
+cc_library(
+    name = "config_bazel_hpp",
+    srcs = ["config/bazel.hpp"],
+    includes = ["config"],
+)
+
+### MFEM Examples #############################################################
+mfem_serial_examples()
+mfem_parallel_examples()
+
+### MFEM library ##############################################################
+
+cc_library(
+    name = "mfem",
+    deps = [
+        "fem",
+        "general",
+        "linalg",
+        "mesh",
+        "@config",
+    ] + select({
+        ":parallel_mode": ["@mpi"],
+        "//conditions:default": [],
+    }),
+)
+
+### Sources ###################################################################
+
+cc_library(
+    name = "fem",
+    srcs = glob([
+        "fem/*.cpp",
+        "fem/ceed/**/*.cpp",
+        "fem/fe/*.cpp",
+        "fem/integ/*.cpp",
+        "fem/lor/*.cpp",
+        # skip moonolith
+        "fem/qinterp/*.cpp",
+        "fem/tmop/*.cpp",
+    ]),
+    deps = [
+        "config_hpp",
+        "fem_hpp",
+        "general_hpp",
+        "linalg_hpp",
+        "mesh_hpp",
+    ] + select({
+        ":parallel_mode": ["@mpi"],
+        "//conditions:default": [],
+    }),
+)
+
+cc_library(
+    name = "general",
+    srcs = glob(["general/*.cpp"]),
+    deps = [
+        "config_hpp",
+        "general_hpp",
+        "linalg_hpp",
+    ] + select({
+        ":parallel_mode": ["@mpi"],
+        "//conditions:default": [],
+    }),
+)
+
+cc_library(
+    name = "linalg",
+    srcs = glob(["linalg/**/*.cpp"]),
+    deps = [
+        "config_hpp",
+        "fem_hpp",
+        "general_hpp",
+        "linalg_hpp",
+        "mesh_hpp",
+    ] + select({
+        ":parallel_mode": ["@mpi"],
+        "//conditions:default": [],
+    }),
+)
+
+cc_library(
+    name = "mesh",
+    srcs = glob([
+        "mesh/*.cpp",
+        "mesh/submesh/*.cpp",
+    ]),
+    deps = [
+        "config_hpp",
+        "fem_hpp",
+        "general_hpp",
+        "linalg_hpp",
+        "mesh_hpp",
+    ] + select({
+        ":parallel_mode": ["@mpi"],
+        "//conditions:default": [],
+    }),
+)
+
+### Headers ###################################################################
+
+cc_library(
+    name = "examples_hpp",
+    srcs = glob(["examples/*.hpp"]),
+)
+
+cc_library(
+    name = "general_hpp",
+    srcs = glob(
+        [
+            "general/*.hpp",
+            "general/*.h",
+        ],
+    ),
+)
+
+cc_library(
+    name = "mfem_hpp",
+    srcs = ["mfem.hpp"],
+    deps = ["config_hpp"],
+)
+
+cc_library(
+    name = "config_hpp",
+    srcs = glob(["config/*.hpp"]),
+    deps = [
+        ":config_bazel_hpp",
+        "@config",
+    ],
+)
+
+cc_library(
+    name = "fem_hpp",
+    srcs = glob([
+        "fem/*.hpp",
+        "fem/*.h",
+    ]),
+    deps = [
+        "config_hpp",
+        "fem_ceed_hpp",
+        "fem_fe_hpp",
+        "fem_integ_hpp",
+        "fem_lor_hpp",
+        # "fem_moonolith_hpp",
+        "fem_qinterp_hpp",
+        "fem_tmop_hpp",
+    ],
+)
+
+cc_library(
+    name = "fem_ceed_hpp",
+    srcs = glob(["fem/ceed/**/*.hpp"]),
+)
+
+cc_library(
+    name = "fem_moonolith_hpp",
+    srcs = glob(["fem/moonolith/**/*.hpp"]),
+)
+
+cc_library(
+    name = "fem_fe_hpp",
+    srcs = glob(["fem/fe/*.hpp"]),
+)
+
+cc_library(
+    name = "fem_integ_hpp",
+    srcs = glob(["fem/integ/*.hpp"]),
+)
+
+cc_library(
+    name = "fem_lor_hpp",
+    srcs = glob(["fem/lor/*.hpp"]),
+)
+
+cc_library(
+    name = "fem_qinterp_hpp",
+    srcs = glob(["fem/qinterp/*.hpp"]),
+)
+
+cc_library(
+    name = "fem_tmop_hpp",
+    srcs = glob(["fem/tmop/*.hpp"]),
+)
+
+cc_library(
+    name = "linalg_hpp",
+    srcs = glob(["linalg/*.hpp"]),
+    deps = [
+        "linalg_batched_hpp",
+        "linalg_simd_hpp",
+    ],
+)
+
+cc_library(
+    name = "linalg_batched_hpp",
+    srcs = glob(["linalg/batched/*.hpp"]),
+)
+
+cc_library(
+    name = "linalg_simd_hpp",
+    srcs = glob(["linalg/simd/*.hpp"]),
+)
+
+cc_library(
+    name = "mesh_hpp",
+    srcs = glob(["mesh/*.hpp"]),
+    deps = ["submesh_hpp"],
+)
+
+cc_library(
+    name = "submesh_hpp",
+    srcs = glob(["mesh/submesh/*.hpp"]),
+)

MODULE.bazel

@@ -0,0 +1,36 @@
+"""
+Bazel Central Registry: https://registry.bazel.build
+"""
+module(name = "mfem", version = "4.7")
+
+# https://github.com/bazelbuild/rules_cc
+# https://registry.bazel.build/modules/rules_cc
+bazel_dep(name = "rules_cc", version = "0.1.1")
+
+# https://github.com/bazel-contrib/toolchains_llvm/releases
+# https://github.com/bazel-contrib/toolchains_llvm/blob/master/toolchain/internal/llvm_distributions.bzl
+bazel_dep(name = "toolchains_llvm", version = "1.3.0")
+
+# Configure and register the toolchain.
+llvm = use_extension("@toolchains_llvm//toolchain/extensions:llvm.bzl", "llvm")
+# llvm.toolchain(llvm_version = "16.0.0")
+llvm.toolchain(llvm_version = "18.1.8")
+# x86_64-linux-gnu-ubuntu-20.04     13.0.0
+# x86_64-linux-gnu-ubuntu-22.04     17.0.6
+# macOS                             17.0.6
+# llvm.toolchain(llvm_version = "17.0.6", stdlib = {"linux-x86_64": "stdc++"}) # libc++ / stdc++
+
+use_repo(llvm, "llvm_toolchain")
+register_toolchains("@llvm_toolchain//:all")
+
+# https://github.com/bazelbuild/bazel-skylib/releases
+bazel_dep(name = "bazel_skylib", version = "1.7.1")
+
+# External dependencies
+external = use_extension("//config/bazel:external.bzl", "external")
+use_repo(external, "mpi")
+use_repo(external, "hypre")
+use_repo(external, "metis")
+use_repo(external, "config")
+
+

config/bazel/config.bzl

@@ -0,0 +1,49 @@
+"""
+Bazel configuration helper functions
+"""
+
+load("@rules_cc//cc:defs.bzl", "cc_binary")
+
+### Examples ##################################################################
+def mfem_serial_examples():
+    for n in range(1, 40):
+        example = "ex" + str(n)
+        cc_binary(
+            name = example,
+            srcs = ["examples/" + example + ".cpp"],
+            deps = [
+                "mfem",
+                "mfem_hpp",
+                "examples_hpp",
+            ],
+        )
+
+def mfem_parallel_examples():
+    for n in range(1, 40):
+        example = "ex" + str(n) + "p"
+        cc_binary(
+            name = example,
+            srcs = ["examples/" + example + ".cpp"],
+            deps = [
+                "mfem",
+                "mfem_hpp",
+                "examples_hpp",
+                "@mpi",
+            ],
+        )
+
+### MFEM_USE_* ################################################################
+def _mfem_use(ctx):
+    value = "//" if not ctx.attr.use else ""
+    value += "#define " + ctx.attr.define
+    return [
+        platform_common.TemplateVariableInfo({ctx.attr.define: value}),
+    ]
+
+mfem_use = rule(
+    implementation = _mfem_use,
+    attrs = {
+        "define": attr.string(),
+        "use": attr.bool(),
+    },
+)

config/bazel/external.bzl

@@ -0,0 +1,57 @@
+"""Modules for dependencies not included in the Bazel Central Registry"""
+# https://bazel.build/rules/lib/repo/local
+
+load("@bazel_tools//tools/build_defs/repo:local.bzl", "new_local_repository")
+
+def _external(_):
+    # MFEM configuration file
+    new_local_repository(
+        name = "config",
+        build_file_content = """
+cc_library(
+    name = "config",
+    defines = ["MFEM_CONFIG_FILE=\\\\\\"config/bazel.hpp\\\\\\""],
+    visibility = ["//visibility:public"],
+)""",
+        path = "",
+    )
+    # MPI implementation
+    new_local_repository(
+        name = "mpi",
+        build_file_content = """
+cc_library(
+    name = "mpi",
+    srcs = ["lib/libmpi.dylib"],
+    hdrs = glob(["include/**/*.h"]),
+    includes = ["include"],
+    visibility = ["//visibility:public"],
+    deps = ["@hypre", "@metis"],
+)""",
+        path = "/opt/homebrew/opt/open-mpi",
+    )
+    new_local_repository(
+        name = "hypre",
+        build_file_content = """
+cc_library(
+    name = "hypre",
+    srcs = ["lib/libHYPRE.a"],
+    hdrs = glob(["include/*.h"]),
+    includes = ["include"],
+    visibility = ["//visibility:public"],
+)""",
+        path = "/opt/homebrew/opt/hypre",
+    )
+    new_local_repository(
+        name = "metis",
+        build_file_content = """
+cc_library(
+    name = "metis",
+    srcs = ["lib/libmetis.dylib"],
+    hdrs = glob(["include/metis.h"]),
+    includes = ["include"],
+    visibility = ["//visibility:public"],
+)""",
+        path = "/opt/homebrew/opt/metis",
+    )
+
+external = module_extension(implementation = _external)

config/bazel/external_linux.bzl

@@ -0,0 +1,61 @@
+"""Modules for dependencies not included in the Bazel Central Registry"""
+# https://bazel.build/rules/lib/repo/local
+
+load("@bazel_tools//tools/build_defs/repo:local.bzl", "new_local_repository")
+
+def _external(_):
+    new_local_repository(
+        name = "config",
+        build_file_content = """
+cc_library(
+    name = "config",
+    defines = ["MFEM_CONFIG_FILE=\\\\\\"config/bazel.hpp\\\\\\""],
+    visibility = ["//visibility:public"],
+)""",
+        path = "",
+    )
+    new_local_repository(
+        name = "mpi",
+        build_file_content = """
+cc_library(
+    name = "mpi",
+    # srcs = ["lib/libmpi.dylib"],
+    srcs = ["lib/libmpi.so"],
+    hdrs = glob(["include/**/*.h"]),
+    includes = ["include"],
+    visibility = ["//visibility:public"],
+    deps = ["@hypre", "@metis"],
+)""",
+        path = "/usr/lib/x86_64-linux-gnu/openmpi"
+    )
+    new_local_repository(
+        name = "hypre",
+        build_file_content = """
+cc_library(
+    name = "hypre",
+    # srcs = ["lib/libHYPRE.a"],
+    srcs = ["lib/x86_64-linux-gnu/libHYPRE.so",
+            "lib/x86_64-linux-gnu/libHYPRE_core.so"],
+    hdrs = glob(["include/hypre/*.h"]),
+    # includes = ["include"],
+    includes = ["include/hypre"],
+    visibility = ["//visibility:public"],
+)""",
+        path = "/usr"
+    )
+    new_local_repository(
+        name = "metis",
+        build_file_content = """
+cc_library(
+    name = "metis",
+    # srcs = ["lib/libmetis.dylib"],
+    srcs = ["lib/x86_64-linux-gnu/libmetis.so.5"],
+    hdrs = glob(["include/metis.h"]),
+    # includes = ["include"],
+    visibility = ["//visibility:public"],
+)""",
+        # path = "/opt/homebrew/Cellar/metis/5.1.0",
+        path = "/usr",
+    )
+
+external = module_extension(implementation = _external)

config/bazel/external_macos.bzl

@@ -0,0 +1,55 @@
+"""Modules for dependencies not included in the Bazel Central Registry"""
+# https://bazel.build/rules/lib/repo/local
+
+load("@bazel_tools//tools/build_defs/repo:local.bzl", "new_local_repository")
+
+def _external(_):
+    new_local_repository(
+        name = "config",
+        build_file_content = """
+cc_library(
+    name = "config",
+    defines = ["MFEM_CONFIG_FILE=\\\\\\"config/bazel.hpp\\\\\\""],
+    visibility = ["//visibility:public"],
+)""",
+        path = "",
+    )
+    new_local_repository(
+        name = "mpi",
+        build_file_content = """
+cc_library(
+    name = "mpi",
+    srcs = ["lib/libmpi.dylib"],
+    hdrs = glob(["include/**/*.h"]),
+    includes = ["include"],
+    visibility = ["//visibility:public"],
+    deps = ["@hypre", "@metis"],
+)""",
+        path = "/opt/homebrew/opt/open-mpi",
+    )
+    new_local_repository(
+        name = "hypre",
+        build_file_content = """
+cc_library(
+    name = "hypre",
+    srcs = ["lib/libHYPRE.a"],
+    hdrs = glob(["include/*.h"]),
+    includes = ["include"],
+    visibility = ["//visibility:public"],
+)""",
+        path = "/opt/homebrew/opt/hypre",
+    )
+    new_local_repository(
+        name = "metis",
+        build_file_content = """
+cc_library(
+    name = "metis",
+    srcs = ["lib/libmetis.dylib"],
+    hdrs = glob(["include/metis.h"]),
+    includes = ["include"],
+    visibility = ["//visibility:public"],
+)""",
+        path = "/opt/homebrew/opt/metis",
+    )
+
+external = module_extension(implementation = _external)

config/bazel/settings.bzl

@@ -0,0 +1,21 @@
+"""
+This module defines build settings: 'mode' and 'precision'.
+"""
+
+# Mode ########################################################################
+ModeInfo = provider(doc = "serial or parallel", fields = ["type"])
+def PrintMode(ctx):
+    ctx.actions.write(output = ctx.outputs.log, 
+                      content = "Compiling in " + ctx.attr.mode[ModeInfo].type + "!")
+print_mode = rule(implementation = PrintMode, attrs = {"mode": attr.label()})
+def Mode(ctx): return ModeInfo(type = ctx.label.name)
+mode = rule(implementation = Mode)
+
+# Precision ###################################################################
+PrecisionInfo = provider(doc = "single or double", fields = ["type"])
+def PrintPrecision(ctx):
+    ctx.actions.write(output = ctx.outputs.log, 
+                      content = "Compiling in " + ctx.attr.precision[PrecisionInfo].type + "!")
+print_precision = rule(implementation = PrintPrecision, attrs = {"precision": attr.label()})
+def Precision(ctx): return PrecisionInfo(type = ctx.label.name)
+precision = rule(implementation = Precision)

fem/bilininteg.cpp

@@ -2829,17 +2829,18 @@ void VectorDivergenceIntegrator::AssembleElementMatrix2(
    DenseMatrix &elmat)
 {
    dim  = trial_fe.GetDim();
+   sdim = Trans.GetSpaceDim();
    int trial_dof = trial_fe.GetDof();
    int test_dof = test_fe.GetDof();
    real_t c;
 
    dshape.SetSize (trial_dof, dim);
-   gshape.SetSize (trial_dof, dim);
-   Jadj.SetSize (dim);
-   divshape.SetSize (dim*trial_dof);
+   gshape.SetSize (trial_dof, sdim);
+   Jadj.SetSize (dim, sdim);
+   divshape.SetSize (sdim*trial_dof);
    shape.SetSize (test_dof);
 
-   elmat.SetSize (test_dof, dim*trial_dof);
+   elmat.SetSize (test_dof, sdim*trial_dof);
 
    const IntegrationRule *ir = GetIntegrationRule(trial_fe, test_fe, Trans);
 
@@ -2853,13 +2854,15 @@ void VectorDivergenceIntegrator::AssembleElementMatrix2(
       trial_fe.CalcDShape (ip, dshape);
       test_fe.CalcPhysShape (Trans, shape);
 
+      // AdjugateJacobian = / adj(J),         if J is square
+      //                    \ adj(J^t.J).J^t, otherwise
       CalcAdjugate(Trans.Jacobian(), Jadj);
-
       Mult (dshape, Jadj, gshape);
 
       gshape.GradToDiv (divshape);
 
       c = ip.weight;
+      if (dim != sdim) { c /= Trans.Weight(); }
       if (Q)
       {
          c *= Q -> Eval (Trans, ip);

fem/bilininteg.hpp

@@ -2935,7 +2935,7 @@ private:
    Vector pa_data;
    const DofToQuad *trial_maps, *test_maps; ///< Not owned
    const GeometricFactors *geom;            ///< Not owned
-   int dim, ne, nq;
+   int dim, sdim, ne, nq;
    int trial_dofs1D, test_dofs1D, quad1D;
 
 public:

fem/tmop.cpp

@@ -3899,8 +3899,18 @@ void TMOP_Integrator::ParUpdateAfterMeshTopologyChange()
 
 real_t TMOP_Integrator::GetElementEnergy(const FiniteElement &el,
                                          ElementTransformation &T,
-                                         const Vector &elfun)
+                                         const Vector &d_el)
 {
+   // Form the Vector of node positions, depending on what's the input.
+   Vector elfun;
+   if (x_0)
+   {
+      // The input is the displacement.
+      x_0->GetElementDofValues(T.ElementNo, elfun);
+      elfun += d_el;
+   }
+   else { elfun = d_el; }
+
    const int dof = el.GetDof(), dim = el.GetDim();
    const int el_id = T.ElementNo;
    real_t energy;
@@ -4197,38 +4207,48 @@ real_t TMOP_Integrator::GetDerefinementElementEnergy(const FiniteElement &el,
 
 void TMOP_Integrator::AssembleElementVector(const FiniteElement &el,
                                             ElementTransformation &T,
-                                            const Vector &elfun, Vector &elvect)
+                                            const Vector &d_el, Vector &elvect)
 {
    if (!fdflag)
    {
-      AssembleElementVectorExact(el, T, elfun, elvect);
+      AssembleElementVectorExact(el, T, d_el, elvect);
    }
    else
    {
-      AssembleElementVectorFD(el, T, elfun, elvect);
+      AssembleElementVectorFD(el, T, d_el, elvect);
    }
 }
 
 void TMOP_Integrator::AssembleElementGrad(const FiniteElement &el,
                                           ElementTransformation &T,
-                                          const Vector &elfun,
+                                          const Vector &d_el,
                                           DenseMatrix &elmat)
 {
    if (!fdflag)
    {
-      AssembleElementGradExact(el, T, elfun, elmat);
+      AssembleElementGradExact(el, T, d_el, elmat);
    }
    else
    {
-      AssembleElementGradFD(el, T, elfun, elmat);
+      AssembleElementGradFD(el, T, d_el, elmat);
    }
 }
 
 void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
                                                  ElementTransformation &T,
-                                                 const Vector &elfun,
+                                                 const Vector &d_el,
                                                  Vector &elvect)
 {
+   // Form the Vector of node positions, depending on what's the input.
+   Vector elfun;
+   if (x_0)
+   {
+      // The input is the displacement.
+      x_0->GetElementDofValues(T.ElementNo, elfun);
+      elfun += d_el;
+   }
+   else { elfun = d_el; }
+
    const int dof = el.GetDof(), dim = el.GetDim();
 
    DenseMatrix Amat(dim), work1(dim), work2(dim);
@@ -4381,9 +4401,19 @@ void TMOP_Integrator::AssembleElementVectorExact(const FiniteElement &el,
 
 void TMOP_Integrator::AssembleElementGradExact(const FiniteElement &el,
                                                ElementTransformation &T,
-                                               const Vector &elfun,
+                                               const Vector &d_el,
                                                DenseMatrix &elmat)
 {
+   // Form the Vector of node positions, depending on what's the input.
+   Vector elfun;
+   if (x_0)
+   {
+      // The input is the displacement.
+      x_0->GetElementDofValues(T.ElementNo, elfun);
+      elfun += d_el;
+   }
+   else { elfun = d_el; }
+
    const int dof = el.GetDof(), dim = el.GetDim();
 
    DSh.SetSize(dof, dim);
@@ -4782,16 +4812,16 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
 
 real_t TMOP_Integrator::GetFDDerivative(const FiniteElement &el,
                                         ElementTransformation &T,
-                                        Vector &elfun, const int dofidx,
+                                        Vector &d_el, const int dofidx,
                                         const int dir, const real_t e_fx,
                                         bool update_stored)
 {
    int dof = el.GetDof();
    int idx = dir*dof+dofidx;
-   elfun[idx]    += dx;
-   real_t e_fxph  = GetElementEnergy(el, T, elfun);
-   elfun[idx]    -= dx;
-   real_t dfdx    = (e_fxph-e_fx)/dx;
+   d_el[idx]     += dx;
+   real_t e_fxph = GetElementEnergy(el, T, d_el);
+   d_el[idx]     -= dx;
+   real_t dfdx   = (e_fxph - e_fx) / dx;
 
    if (update_stored)
    {
@@ -4804,11 +4834,21 @@ real_t TMOP_Integrator::GetFDDerivative(const FiniteElement &el,
 
 void TMOP_Integrator::AssembleElementVectorFD(const FiniteElement &el,
                                               ElementTransformation &T,
-                                              const Vector &elfun,
+                                              const Vector &d_el,
                                               Vector &elvect)
 {
+   // Form the Vector of node positions, depending on what's the input.
+   Vector elfun;
+   if (x_0)
+   {
+      // The input is the displacement.
+      x_0->GetElementDofValues(T.ElementNo, elfun);
+      elfun += d_el;
+   }
+   else { elfun = d_el; }
+
    const int dof = el.GetDof(), dim = el.GetDim(), elnum = T.ElementNo;
-   if (elnum>=ElemDer.Size())
+   if (elnum >= ElemDer.Size())
    {
       ElemDer.Append(new Vector);
       ElemPertEnergy.Append(new Vector);
@@ -4817,14 +4857,14 @@ void TMOP_Integrator::AssembleElementVectorFD(const FiniteElement &el,
    }
 
    elvect.SetSize(dof*dim);
-   Vector elfunmod(elfun);
 
    // In GetElementEnergy(), skip terms that have exact derivative calculations.
    fd_call_flag = true;
 
    // Energy for unperturbed configuration.
-   const real_t e_fx = GetElementEnergy(el, T, elfun);
+   const real_t e_fx = GetElementEnergy(el, T, d_el);
 
+   Vector d_el_mod(d_el);
    for (int j = 0; j < dim; j++)
    {
       for (int i = 0; i < dof; i++)
@@ -4834,7 +4874,7 @@ void TMOP_Integrator::AssembleElementVectorFD(const FiniteElement &el,
             discr_tc->UpdateTargetSpecificationAtNode(
                el, T, i, j, discr_tc->GetTspecPert1H());
          }
-         elvect(j*dof+i) = GetFDDerivative(el, T, elfunmod, i, j, e_fx, true);
+         elvect(j*dof+i) = GetFDDerivative(el, T, d_el_mod, i, j, e_fx, true);
          if (discr_tc) { discr_tc->RestoreTargetSpecificationAtNode(T, i); }
       }
    }
@@ -4872,18 +4912,28 @@ void TMOP_Integrator::AssembleElementVectorFD(const FiniteElement &el,
 
 void TMOP_Integrator::AssembleElementGradFD(const FiniteElement &el,
                                             ElementTransformation &T,
-                                            const Vector &elfun,
+                                            const Vector &d_el,
                                             DenseMatrix &elmat)
 {
+   // Form the Vector of node positions, depending on what's the input.
+   Vector elfun;
+   if (x_0)
+   {
+      // The input is the displacement.
+      x_0->GetElementDofValues(T.ElementNo, elfun);
+      elfun += d_el;
+   }
+   else { elfun = d_el; }
+
    const int dof = el.GetDof(), dim = el.GetDim();
 
    elmat.SetSize(dof*dim);
-   Vector elfunmod(elfun);
 
    const Vector &ElemDerLoc = *(ElemDer[T.ElementNo]);
    const Vector &ElemPertLoc = *(ElemPertEnergy[T.ElementNo]);
 
    // In GetElementEnergy(), skip terms that have exact derivative calculations.
+   Vector d_el_mod(d_el);
    fd_call_flag = true;
    for (int i = 0; i < dof; i++)
    {
@@ -4893,7 +4943,7 @@ void TMOP_Integrator::AssembleElementGradFD(const FiniteElement &el,
          {
             for (int k2 = 0; k2 < dim; k2++)
             {
-               elfunmod(k2*dof+j) += dx;
+               d_el_mod(k2 * dof + j) += dx;
 
                if (discr_tc)
                {
@@ -4920,10 +4970,10 @@ void TMOP_Integrator::AssembleElementGradFD(const FiniteElement &el,
                   }
                }
 
-               real_t e_fx   = ElemPertLoc(k2*dof+j);
-               real_t e_fpxph = GetFDDerivative(el, T, elfunmod, i, k1, e_fx,
+               real_t e_fx    = ElemPertLoc(k2 * dof + j);
+               real_t e_fpxph = GetFDDerivative(el, T, d_el_mod, i, k1, e_fx,
                                                 false);
-               elfunmod(k2*dof+j) -= dx;
+               d_el_mod(k2 * dof + j) -= dx;
                real_t e_fpx = ElemDerLoc(k1*dof+i);
 
                elmat(k1*dof+i, k2*dof+j) = (e_fpxph - e_fpx) / dx;
@@ -5016,6 +5066,22 @@ void TMOP_Integrator::ParEnableNormalization(const ParGridFunction &x)
 }
 #endif
 
+void TMOP_Integrator::SetInitialMeshPos(const GridFunction *x0)
+{
+   x_0 = x0;
+
+   // Compute PA.X0 when we're setting x_0 to something.
+   // TODO move(or copy?) this in AssemblePA.
+   if (PA.enabled && x_0 != nullptr)
+   {
+      const ElementDofOrdering ord = ElementDofOrdering::LEXICOGRAPHIC;
+      const Operator *n0_R = x0->FESpace()->GetElementRestriction(ord);
+      PA.X0.SetSize(n0_R->Height(), Device::GetMemoryType());
+      PA.X0.UseDevice(true);
+      n0_R->Mult(*x_0, PA.X0);
+   }
+}
+
 void TMOP_Integrator::ComputeNormalizationEnergies(const GridFunction &x,
                                                    real_t &metric_energy,
                                                    real_t &lim_energy,

fem/tmop.hpp

@@ -1887,6 +1887,15 @@ protected:
    friend class TMOPNewtonSolver;
    friend class TMOPComboIntegrator;
 
+   // Initial positions of the mesh nodes. Not owned. The pointer is set at the
+   // start of the solve by TMOPNewtonSolver::Mult(), and unset at the end.
+   // When x_0 == nullptr, the integrator works on the mesh positions.
+   // When x_0 != nullptr, the integrator works on the displacements.
+   // TODO in MFEM-5.0 make it always work with displacements.
+   const GridFunction *x_0;
+   // Called with nullptr to unset the x_0 after the problem is solved.
+   void SetInitialMeshPos(const GridFunction *x0);
+
    TMOP_QualityMetric *h_metric;
    TMOP_QualityMetric *metric;        // not owned
    const TargetConstructor *targetC;  // not owned
@@ -1973,7 +1982,9 @@ protected:
    //  E: Q-vector for TMOP-energy
    //     Used as temporary storage when the total energy is computed.
    //  O: Q-Vector of 1.0, used to compute sums using the dot product kernel.
-   // X0: E-vector for initial nodal coordinates used for limiting.
+   // X0: E-vector for initial nodal coordinates.
+   //     Does not change during the TMOP iteration.
+   // XL: E-vector for nodal coordinates used for limiting.
    //     Does not change during the TMOP iteration.
    //  H: Q-Vector for Hessian associated with the metric term.
    //     Updated by every call to PANonlinearFormExtension::GetGradient().
@@ -2006,7 +2017,7 @@ protected:
       mutable DenseTensor Jtr;
       mutable bool Jtr_needs_update;
       mutable bool Jtr_debug_grad;
-      mutable Vector E, O, X0, H, C0, LD, H0, MC;
+      mutable Vector E, O, X0, XL, H, C0, LD, H0, MC;
       const DofToQuad *maps;
       const DofToQuad *maps_lim = nullptr;
       const GeometricFactors *geom;
@@ -2020,20 +2031,20 @@ protected:
 
    void AssembleElementVectorExact(const FiniteElement &el,
                                    ElementTransformation &T,
-                                   const Vector &elfun, Vector &elvect);
+                                   const Vector &d_el, Vector &elvect);
 
    void AssembleElementGradExact(const FiniteElement &el,
                                  ElementTransformation &T,
-                                 const Vector &elfun, DenseMatrix &elmat);
+                                 const Vector &d_el, DenseMatrix &elmat);
 
    void AssembleElementVectorFD(const FiniteElement &el,
                                 ElementTransformation &T,
-                                const Vector &elfun, Vector &elvect);
+                                const Vector &d_el, Vector &elvect);
 
    // Assumes that AssembleElementVectorFD has been called.
    void AssembleElementGradFD(const FiniteElement &el,
                               ElementTransformation &T,
-                              const Vector &elfun, DenseMatrix &elmat);
+                              const Vector &d_el, DenseMatrix &elmat);
 
    void AssembleElemVecAdaptLim(const FiniteElement &el,
                                 IsoparametricTransformation &Tpr,
@@ -2056,7 +2067,7 @@ protected:
 
    real_t GetFDDerivative(const FiniteElement &el,
                           ElementTransformation &T,
-                          Vector &elfun, const int nodenum,const int idir,
+                          Vector &d_el, const int nodenum, const int idir,
                           const real_t baseenergy, bool update_stored);
 
    /** @brief Determines the perturbation, h, for FD-based approximation. */
@@ -2141,7 +2152,7 @@ public:
        @param[in] hm   TMOP_QualityMetric for h-adaptivity (not owned). */
    TMOP_Integrator(TMOP_QualityMetric *m, TargetConstructor *tc,
                    TMOP_QualityMetric *hm)
-      : h_metric(hm), metric(m), targetC(tc), IntegRules(NULL),
+      : x_0(nullptr), h_metric(hm), metric(m), targetC(tc), IntegRules(NULL),
         integ_order(-1), metric_coeff(NULL), metric_normal(1.0),
         lim_nodes0(NULL), lim_coeff(NULL),
         lim_dist(NULL), lim_func(NULL), lim_normal(1.0),
@@ -2323,10 +2334,10 @@ public:
    /** @brief Computes the integral of W(Jacobian(Trt)) over a target zone.
        @param[in] el     Type of FiniteElement.
        @param[in] T      Mesh element transformation.
-       @param[in] elfun  Physical coordinates of the zone. */
+       @param[in] d_el   Physical displacement of the zone w.r.t. x_0. */
    real_t GetElementEnergy(const FiniteElement &el,
                            ElementTransformation &T,
-                           const Vector &elfun) override;
+                           const Vector &d_el) override;
 
    /** @brief Computes the mean of the energies of the given element's children.
 
@@ -2346,11 +2357,11 @@ public:
 
    void AssembleElementVector(const FiniteElement &el,
                               ElementTransformation &T,
-                              const Vector &elfun, Vector &elvect) override;
+                              const Vector &d_el, Vector &elvect) override;
 
    void AssembleElementGrad(const FiniteElement &el,
                             ElementTransformation &T,
-                            const Vector &elfun, DenseMatrix &elmat) override;
+                            const Vector &d_el, DenseMatrix &elmat) override;
 
    TMOP_QualityMetric &GetAMRQualityMetric() { return *h_metric; }
 
@@ -2363,7 +2374,7 @@ public:
    using NonlinearFormIntegrator::AssemblePA;
    void AssemblePA(const FiniteElementSpace&) override;
 
-   void AssembleGradPA(const Vector&, const FiniteElementSpace&) override;
+   void AssembleGradPA(const Vector &, const FiniteElementSpace &) override;
 
    real_t GetLocalStateEnergyPA(const Vector&) const override;
 
@@ -2411,9 +2422,17 @@ public:
 class TMOPComboIntegrator : public NonlinearFormIntegrator
 {
 protected:
+   friend class TMOPNewtonSolver;
+
    // Integrators in the combination. Owned.
    Array<TMOP_Integrator *> tmopi;
 
+   void SetInitialMeshPos(const GridFunction *x0)
+   {
+      for (int i = 0; i < tmopi.Size(); i++)
+      { tmopi[i]->SetInitialMeshPos(x0); }
+   }
+
 public:
    TMOPComboIntegrator() : tmopi(0) { }
 

fem/tmop/tmop_pa.cpp

@@ -20,7 +20,7 @@
 namespace mfem
 {
 
-void TMOP_Integrator::AssembleGradPA(const Vector &xe,
+void TMOP_Integrator::AssembleGradPA(const Vector &de,
                                      const FiniteElementSpace &fes)
 {
    MFEM_VERIFY(PA.enabled, "PA extension setup has not been done!");
@@ -29,6 +29,15 @@ void TMOP_Integrator::AssembleGradPA(const Vector &xe,
    // AssemblePA() was called has not been modified or completely destroyed and
    // a new object created at the same address.
 
+   // Form the Vector of node positions, depending on what's the input.
+   Vector xe(de.Size());
+   if (x_0)
+   {
+      // The input is the displacement.
+      add(PA.X0, de, xe);
+   }
+   else { xe = de; }
+
    if (PA.Jtr_needs_update || targetC->UsesPhysicalCoordinates())
    {
       ComputeAllElementTargets(xe);
@@ -96,12 +105,12 @@ void TMOP_Integrator::AssemblePA_Limiting()
       }
    }
 
-   // lim_nodes0 -> PA.X0 (E-vector)
+   // lim_nodes0 -> PA.XL (E-vector)
    MFEM_VERIFY(lim_nodes0->FESpace() == fes, "");
    const Operator *n0_R = fes->GetElementRestriction(ordering);
-   PA.X0.SetSize(n0_R->Height(), Device::GetMemoryType());
-   PA.X0.UseDevice(true);
-   n0_R->Mult(*lim_nodes0, PA.X0);
+   PA.XL.SetSize(n0_R->Height(), Device::GetMemoryType());
+   PA.XL.UseDevice(true);
+   n0_R->Mult(*lim_nodes0, PA.XL);
 
    // Limiting distances: lim_dist -> PA.LD (E-vector)
    // TODO: remove the hack for the case lim_dist == NULL.
@@ -115,10 +124,7 @@ void TMOP_Integrator::AssemblePA_Limiting()
       const Operator *ld_R = limfes->GetElementRestriction(ordering);
       ld_R->Mult(*lim_dist, PA.LD);
    }
-   else
-   {
-      PA.LD = 1.0;
-   }
+   else { PA.LD = 1.0; }
 }
 
 void TargetConstructor::ComputeAllElementTargets(const FiniteElementSpace &fes,
@@ -221,19 +227,24 @@ void TMOP_Integrator::AssemblePA(const FiniteElementSpace &fes)
    Mesh *mesh = fes.GetMesh();
    const int ne = PA.ne = mesh->GetNE();
    const int dim = PA.dim = mesh->Dimension();
+
    MFEM_VERIFY(PA.dim == 2 || PA.dim == 3, "Not yet implemented!");
    MFEM_VERIFY(mesh->GetNumGeometries(dim) <= 1,
-               "mixed meshes are not supported");
+               "TMOP+PA does not support mixed meshes.");
+   MFEM_VERIFY(mesh->HasGeometry(Geometry::SQUARE) ||
+               mesh->HasGeometry(Geometry::CUBE),
+               "TMOP+PA only supports squares and cubes.");
    MFEM_VERIFY(!fes.IsVariableOrder(), "variable orders are not supported");
+   MFEM_VERIFY(fes.GetOrdering() == Ordering::byNODES,
+               "TMOP+PAP only supports Ordering::byNODES!");
+
    const FiniteElement &fe = *fes.GetTypicalFE();
    PA.ir = &EnergyIntegrationRule(fe);
    const IntegrationRule &ir = *PA.ir;
-   MFEM_VERIFY(fes.GetOrdering() == Ordering::byNODES,
-               "PA Only supports Ordering::byNODES!");
-
    const int nq = PA.nq = ir.GetNPoints();
    const DofToQuad::Mode mode = DofToQuad::TENSOR;
    PA.maps = &fe.GetDofToQuad(ir, mode);
+   // Note - initial mesh. TODO delete this?
    PA.geom = mesh->GetGeometricFactors(ir, GeometricFactors::JACOBIANS);
 
    // Energy vector, scalar Q-vector
@@ -265,6 +276,7 @@ void TMOP_Integrator::AssemblePA(const FiniteElementSpace &fes)
             ElementTransformation& T = *PA.fes->GetElementTransformation(e);
             for (int q = 0; q < ir.GetNPoints(); ++q)
             {
+               // Note that this is always on the initial mesh.
                M0(q,e) = metric_coeff->Eval(T, ir.IntPoint(q));
             }
          }
@@ -282,7 +294,7 @@ void TMOP_Integrator::AssemblePA(const FiniteElementSpace &fes)
    PA.Jtr_needs_update = true;
    PA.Jtr_debug_grad = false;
 
-   // Limiting: lim_coeff -> PA.C0, lim_nodes0 -> PA.X0, lim_dist -> PA.LD, PA.H0
+   // Limiting: lim_coeff -> PA.C0, lim_nodes0 -> PA.XL, lim_dist -> PA.LD, PA.H0
    if (lim_coeff) { AssemblePA_Limiting(); }
 }
 
@@ -311,10 +323,19 @@ void TMOP_Integrator::AssembleGradDiagonalPA(Vector &de) const
    }
 }
 
-void TMOP_Integrator::AddMultPA(const Vector &xe, Vector &ye) const
+void TMOP_Integrator::AddMultPA(const Vector &de, Vector &ye) const
 {
    // This method must be called after AssemblePA().
 
+   // Form the Vector of node positions, depending on what's the input.
+   Vector xe(de.Size());
+   if (x_0)
+   {
+      // The input is the displacement.
+      add(PA.X0, de, xe);
+   }
+   else { xe = de; }
+
    if (PA.Jtr_needs_update || targetC->UsesPhysicalCoordinates())
    {
       ComputeAllElementTargets(xe);
@@ -358,10 +379,19 @@ void TMOP_Integrator::AddMultGradPA(const Vector &re, Vector &ce) const
    }
 }
 
-real_t TMOP_Integrator::GetLocalStateEnergyPA(const Vector &xe) const
+real_t TMOP_Integrator::GetLocalStateEnergyPA(const Vector &de) const
 {
    // This method must be called after AssemblePA().
 
+   // Form the Vector of node positions, depending on what's the input.
+   Vector xe(de.Size());
+   if (x_0)
+   {
+      // The input is the displacement.
+      add(PA.X0, de, xe);
+   }
+   else { xe = de; }
+
    real_t energy = 0.0;
 
    if (PA.Jtr_needs_update || targetC->UsesPhysicalCoordinates())

fem/tmop/tmop_pa_h2s_c0.cpp

@@ -161,13 +161,13 @@ void TMOP_Integrator::AssembleGradPA_C0_2D(const Vector &X) const
    const Array<real_t> &B   = PA.maps->B;
    const Array<real_t> &BLD = PA.maps_lim->B;
    const Vector &C0 = PA.C0;
-   const Vector &X0 = PA.X0;
+   const Vector &XL = PA.XL;
    Vector &H0 = PA.H0;
 
    auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
    const bool exp_lim = (el) ? true : false;
 
-   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,H0,
+   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,XL,X,H0,
                            exp_lim);
 }
 

fem/tmop/tmop_pa_h3s_c0.cpp

@@ -179,13 +179,13 @@ void TMOP_Integrator::AssembleGradPA_C0_3D(const Vector &X) const
    const Array<real_t> &B   = PA.maps->B;
    const Array<real_t> &BLD = PA.maps_lim->B;
    const Vector &C0 = PA.C0;
-   const Vector &X0 = PA.X0;
+   const Vector &XL = PA.XL;
    Vector &H0 = PA.H0;
 
    auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
    const bool exp_lim = (el) ? true : false;
 
-   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,
+   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,XL,X,
                            H0,exp_lim);
 }
 

fem/tmop/tmop_pa_p2_c0.cpp

@@ -163,12 +163,12 @@ void TMOP_Integrator::AddMultPA_C0_2D(const Vector &X, Vector &Y) const
    const Array<real_t> &BLD = PA.maps_lim->B;
    MFEM_VERIFY(PA.maps_lim->ndof == D1D, "");
    MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
-   const Vector &X0 = PA.X0;
+   const Vector &XL = PA.XL;
    const Vector &C0 = PA.C0;
    auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
    const bool exp_lim = (el) ? true : false;
 
-   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y,
+   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,XL,X,Y,
                            exp_lim);
 }
 

fem/tmop/tmop_pa_p3_c0.cpp

@@ -171,12 +171,12 @@ void TMOP_Integrator::AddMultPA_C0_3D(const Vector &X, Vector &Y) const
    const Array<real_t> &BLD = PA.maps_lim->B;
    MFEM_VERIFY(PA.maps_lim->ndof == D1D, "");
    MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
-   const Vector &X0 = PA.X0;
+   const Vector &XL = PA.XL;
    const Vector &C0 = PA.C0;
    auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
    const bool exp_lim = (el) ? true : false;
 
-   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y,
+   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,XL,X,Y,
                            exp_lim);
 }
 

fem/tmop/tmop_pa_w2_c0.cpp

@@ -142,7 +142,7 @@ real_t TMOP_Integrator::GetLocalStateEnergyPA_C0_2D(const Vector &X) const
    const Array<real_t> &BLD = PA.maps_lim->B;
    MFEM_VERIFY(PA.maps_lim->ndof == D1D, "");
    MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
-   const Vector &X0 = PA.X0;
+   const Vector &XL = PA.XL;
    const Vector &C0 = PA.C0;
    const Vector &O = PA.O;
    Vector &E = PA.E;
@@ -150,7 +150,7 @@ real_t TMOP_Integrator::GetLocalStateEnergyPA_C0_2D(const Vector &X) const
    auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
    const bool exp_lim = (el) ? true : false;
 
-   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E,
+   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,XL,X,O,E,
                            exp_lim);
 }
 

fem/tmop/tmop_pa_w3_c0.cpp

@@ -154,7 +154,7 @@ real_t TMOP_Integrator::GetLocalStateEnergyPA_C0_3D(const Vector &X) const
    const Array<real_t> &BLD = PA.maps_lim->B;
    MFEM_VERIFY(PA.maps_lim->ndof == D1D, "");
    MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
-   const Vector &X0 = PA.X0;
+   const Vector &XL = PA.XL;
    const Vector &C0 = PA.C0;
    const Vector &O = PA.O;
    Vector &E = PA.E;
@@ -162,7 +162,7 @@ real_t TMOP_Integrator::GetLocalStateEnergyPA_C0_3D(const Vector &X) const
    auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
    const bool exp_lim = (el) ? true : false;
 
-   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E,
+   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,XL,X,O,E,
                            exp_lim);
 }
 

fem/tmop_tools.cpp

@@ -417,9 +417,12 @@ void InterpolatorFP::ComputeAtGivenPositions(const Vector &positions,
 
 #endif
 
-real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
+real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &d_in,
                                               const Vector &b) const
 {
+   Vector x_in(x_0.Size());
+   add(x_0, d_in, x_in);
+
    const FiniteElementSpace *fes = NULL;
    real_t energy_in = 0.0;
 #ifdef MFEM_USE_MPI
@@ -428,7 +431,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
    if (parallel)
    {
       fes = p_nlf->FESpace();
-      energy_in = p_nlf->GetEnergy(x);
+      energy_in = p_nlf->GetEnergy(d_in);
    }
 #endif
    const bool serial = !parallel;
@@ -437,7 +440,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
    if (serial)
    {
       fes = nlf->FESpace();
-      energy_in = nlf->GetEnergy(x);
+      energy_in = nlf->GetEnergy(d_in);
    }
 
    // Get the local prolongation of the solution vector.
@@ -446,13 +449,13 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
    if (serial)
    {
       const SparseMatrix *cP = fes->GetConformingProlongation();
-      if (!cP) { x_out_loc = x; }
-      else     { cP->Mult(x, x_out_loc); }
+      if (!cP) { x_out_loc = x_in; }
+      else     { cP->Mult(x_in, x_out_loc); }
    }
 #ifdef MFEM_USE_MPI
    else
    {
-      fes->GetProlongationMatrix()->Mult(x, x_out_loc);
+      fes->GetProlongationMatrix()->Mult(x_in, x_out_loc);
    }
 #endif
 
@@ -505,7 +508,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
 
    const bool have_b = (b.Size() == Height());
 
-   Vector x_out(x.Size());
+   Vector x_out(x_in.Size()), d_out(d_in.Size());
    bool x_out_ok = false;
    real_t energy_out = 0.0, min_detT_out;
    const real_t norm_in = Norm(r);
@@ -523,8 +526,13 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       avg_fit_err = 0.0;
       max_fit_err = 0.0;
 
+      //
       // Update the mesh and get the L-vector in x_out_loc.
-      add(x, -scale, c, x_out);
+      //
+      // Form limited (line-search) displacement d_out = d_in - scale * c,
+      // and the corresponding mesh positions x_out = x_0 + d_out.
+      add(d_in, -scale, c, d_out);
+      add(x_0, d_out, x_out);
       if (serial)
       {
          const SparseMatrix *cP = fes->GetConformingProlongation();
@@ -561,8 +569,8 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       // energy and residual, so their increase/decrease is not relevant.
       if (untangling) { x_out_ok = true; break; }
 
-      // Check the changes in total energy.
-      ProcessNewState(x_out);
+      // Update mesh-dependent quantities.
+      ProcessNewState(d_out);
 
       // Ensure sufficient decrease in fitting error if we are trying to
       // converge based on error.
@@ -579,14 +587,15 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
          }
       }
 
+      // Check the changes in total energy.
       if (serial)
       {
-         energy_out = nlf->GetGridFunctionEnergy(x_out_loc);
+         energy_out = nlf->GetEnergy(d_out);
       }
 #ifdef MFEM_USE_MPI
       else
       {
-         energy_out = p_nlf->GetParGridFunctionEnergy(x_out_loc);
+         energy_out = p_nlf->GetEnergy(d_out);
       }
 #endif
       if (energy_out > energy_in + 0.2*fabs(energy_in) ||
@@ -601,7 +610,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
       }
 
       // Check the changes in the Newton residual.
-      oper->Mult(x_out, r);
+      oper->Mult(d_out, r);
       if (have_b) { r -= b; }
       real_t norm_out = Norm(r);
 
@@ -656,6 +665,49 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
    return scale;
 }
 
+void TMOPNewtonSolver::Mult(const Vector &b, Vector &x) const
+{
+   x_0 = x;
+
+   //
+   // Pass down the initial position to the integrators.
+   //
+   // Prolongate x to ldofs.
+   const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
+   GridFunction x_0_loc(const_cast<FiniteElementSpace *>(nlf->FESpace()));
+   const Operator *P = nlf->GetProlongation();
+   // TODO if (periodic) { x_0_loc = x }
+   if (P) { P->Mult(x, x_0_loc); }
+   else   { x_0_loc = x; }
+   // Pass the positions to the integrators.
+   const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
+   for (int i = 0; i < integs.Size(); i++)
+   {
+      auto ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+      if (ti) { ti->SetInitialMeshPos(&x_0_loc); }
+      auto co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+      if (co) { co->SetInitialMeshPos(&x_0_loc); }
+   }
+
+   // We solve for the displacement, which always starts from zero.
+   Vector d(x.Size()); d = 0.0;
+   if (solver_type == 0)      { NewtonSolver::Mult(b, d); }
+   else if (solver_type == 1) { LBFGSSolver::Mult(b, d); }
+   else { MFEM_ABORT("Invalid solver_type"); }
+
+   // Form the final mesh using the computed displacement.
+   x += d;
+
+   // Make sure the pointers don't use invalid memory (x_0_loc is gone).
+   for (int i = 0; i < integs.Size(); i++)
+   {
+      auto ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+      if (ti) { ti->SetInitialMeshPos(nullptr); }
+      auto co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+      if (co) { co->SetInitialMeshPos(nullptr); }
+   }
+}
+
 void TMOPNewtonSolver::UpdateSurfaceFittingWeight(real_t factor) const
 {
    const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
@@ -788,8 +840,11 @@ bool TMOPNewtonSolver::IsSurfaceFittingEnabled() const
    return false;
 }
 
-void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
+void TMOPNewtonSolver::ProcessNewState(const Vector &d) const
 {
+   Vector x(x_0.Size());
+   add(x_0, d, x);
+
    const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
    const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
 
@@ -819,30 +874,15 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
    }
 
    Vector x_loc;
-   const FiniteElementSpace *x_fes = nullptr;
-   if (parallel)
+   const Operator *P = nlf->GetProlongation();
+   if (P)
    {
-#ifdef MFEM_USE_MPI
-      const ParNonlinearForm *pnlf =
-         dynamic_cast<const ParNonlinearForm *>(oper);
-
-      x_fes = pnlf->ParFESpace();
-      x_loc.SetSize(x_fes->GetVSize());
-      x_fes->GetProlongationMatrix()->Mult(x, x_loc);
-#endif
-   }
-   else
-   {
-      x_fes = nlf->FESpace();
-      const Operator *P = nlf->GetProlongation();
-      if (P)
-      {
-         x_loc.SetSize(P->Height());
-         P->Mult(x,x_loc);
-      }
-      else { x_loc = x; }
+      x_loc.SetSize(P->Height());
+      P->Mult(x, x_loc);
    }
+   else { x_loc = x; }
 
+   const FiniteElementSpace *x_fes = nlf->FESpace();
    for (int i = 0; i < integs.Size(); i++)
    {
       ti = dynamic_cast<TMOP_Integrator *>(integs[i]);

fem/tmop_tools.hpp

@@ -170,6 +170,10 @@ protected:
    int solver_type;
    bool parallel;
 
+   // Starting mesh positions (tdofs). Updated by the call to Mult().
+   // This solver solves for d, where the final mesh is x = x_0 + d.
+   mutable Vector x_0;
+
    // Line search step is rejected if min(detJ) <= min_detJ_limit.
    real_t min_detJ_limit = 0.0;
 
@@ -235,11 +239,11 @@ protected:
 public:
 #ifdef MFEM_USE_MPI
    TMOPNewtonSolver(MPI_Comm comm, const IntegrationRule &irule, int type = 0)
-      : LBFGSSolver(comm), solver_type(type), parallel(true),
+      : LBFGSSolver(comm), solver_type(type), parallel(true), x_0(0),
         ir(irule), IntegRules(NULL), integ_order(-1) { }
 #endif
    TMOPNewtonSolver(const IntegrationRule &irule, int type = 0)
-      : LBFGSSolver(), solver_type(type), parallel(false),
+      : LBFGSSolver(), solver_type(type), parallel(false), x_0(0),
         ir(irule), IntegRules(NULL), integ_order(-1) { }
 
    /// Prescribe a set of integration rules; relevant for mixed meshes.
@@ -259,11 +263,11 @@ public:
    /// Compute scaling factor for the node movement direction using line-search.
    /// We impose constraints on TMOP energy, gradient, minimum Jacobian of
    /// the mesh, and (optionally) on the surface fitting error.
-   real_t ComputeScalingFactor(const Vector &x, const Vector &b) const override;
+   real_t ComputeScalingFactor(const Vector &d, const Vector &b) const override;
 
    /// Update (i) discrete functions at new nodal positions, and
    /// (ii) surface fitting weight.
-   void ProcessNewState(const Vector &x) const override;
+   void ProcessNewState(const Vector &d) const override;
 
    /** @name Methods for adaptive surface fitting.
        \brief These methods control the behavior of the weight and the
@@ -355,18 +359,8 @@ public:
       min_detJ_limit = threshold;
    }
 
-   void Mult(const Vector &b, Vector &x) const override
-   {
-      if (solver_type == 0)
-      {
-         NewtonSolver::Mult(b, x);
-      }
-      else if (solver_type == 1)
-      {
-         LBFGSSolver::Mult(b, x);
-      }
-      else { MFEM_ABORT("Invalid type"); }
-   }
+   /// Optimizes the mesh positions given by @a x.
+   void Mult(const Vector &b, Vector &x) const override;
 
    void SetSolver(Solver &solver) override
    {

miniapps/meshing/mesh-optimizer.cpp

@@ -152,7 +152,7 @@ int main(int argc, char *argv[])
    int barrier_type       = 0;
    int worst_case_type    = 0;
 
-   // 1. Parse command-line options.
+   // Parse command-line options.
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
                   "Mesh file to use.");
@@ -325,17 +325,17 @@ int main(int argc, char *argv[])
    Device device(devopt);
    device.Print();
 
-   // 2. Initialize and refine the starting mesh.
+   // Initialize and refine the starting mesh.
    Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
    for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
    const int dim = mesh->Dimension();
 
    if (hradaptivity) { mesh->EnsureNCMesh(); }
 
-   // 3. Define a finite element space on the mesh-> Here we use vector finite
-   //    elements which are tensor products of quadratic finite elements. The
-   //    number of components in the vector finite element space is specified by
-   //    the last parameter of the FiniteElementSpace constructor.
+   // Define a finite element space on the mesh-> Here we use vector finite
+   // elements which are tensor products of quadratic finite elements. The
+   // number of components in the vector finite element space is specified by
+   // the last parameter of the FiniteElementSpace constructor.
    FiniteElementCollection *fec;
    if (mesh_poly_deg <= 0)
    {
@@ -346,25 +346,25 @@ int main(int argc, char *argv[])
    FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec, dim,
                                                         mesh_node_ordering);
 
-   // 4. Make the mesh curved based on the above finite element space. This
-   //    means that we define the mesh elements through a fespace-based
-   //    transformation of the reference element.
+   // Make the mesh curved based on the above finite element space. This
+   // means that we define the mesh elements through a fespace-based
+   // transformation of the reference element.
    mesh->SetNodalFESpace(fespace);
 
-   // 5. Set up an empty right-hand side vector b, which is equivalent to b=0.
+   // Set up an empty right-hand side vector b, which is equivalent to b=0.
    Vector b(0);
 
-   // 6. Get the mesh nodes (vertices and other degrees of freedom in the finite
-   //    element space) as a finite element grid function in fespace. Note that
-   //    changing x automatically changes the shapes of the mesh elements.
+   // Get the mesh nodes (vertices and other degrees of freedom in the finite
+   // element space) as a finite element grid function in fespace. Note that
+   // changing x automatically changes the shapes of the mesh elements.
    GridFunction x(fespace);
    mesh->SetNodalGridFunction(&x);
 
-   // 7. Define a vector representing the minimal local mesh size in the mesh
-   //    nodes. We index the nodes using the scalar version of the degrees of
-   //    freedom in fespace. Note: this is partition-dependent.
+   // Define a vector representing the minimal local mesh size in the mesh
+   // nodes. We index the nodes using the scalar version of the degrees of
+   // freedom in fespace. Note: this is partition-dependent.
    //
-   //    In addition, compute average mesh size and total volume.
+   // In addition, compute average mesh size and total volume.
    Vector h0(fespace->GetNDofs());
    h0 = infinity();
    real_t mesh_volume = 0.0;
@@ -383,48 +383,50 @@ int main(int argc, char *argv[])
    }
    const real_t small_phys_size = pow(mesh_volume, 1.0 / dim) / 100.0;
 
-   // 8. Add a random perturbation to the nodes in the interior of the domain.
-   //    We define a random grid function of fespace and make sure that it is
-   //    zero on the boundary and its values are locally of the order of h0.
-   //    The latter is based on the DofToVDof() method which maps the scalar to
-   //    the vector degrees of freedom in fespace.
-   GridFunction rdm(fespace);
-   rdm.Randomize();
-   rdm -= 0.25; // Shift to random values in [-0.5,0.5].
-   rdm *= jitter;
-   rdm.HostReadWrite();
-   // Scale the random values to be of order of the local mesh size.
-   for (int i = 0; i < fespace->GetNDofs(); i++)
+   // Add a random perturbation to the nodes in the interior of the domain.
+   // We define a random grid function of fespace and make sure that it is
+   // zero on the boundary and its values are locally of the order of h0.
+   // The latter is based on the DofToVDof() method which maps the scalar to
+   // the vector degrees of freedom in fespace.
+   if (jitter > 0)
    {
-      for (int d = 0; d < dim; d++)
+      GridFunction rdm(fespace);
+      rdm.Randomize();
+      rdm -= 0.25; // Shift to random values in [-0.5,0.5].
+      rdm *= jitter;
+      rdm.HostReadWrite();
+      // Scale the random values to be of order of the local mesh size.
+      for (int i = 0; i < fespace->GetNDofs(); i++)
       {
-         rdm(fespace->DofToVDof(i,d)) *= h0(i);
+         for (int d = 0; d < dim; d++)
+         {
+            rdm(fespace->DofToVDof(i,d)) *= h0(i);
+         }
       }
+      Array<int> vdofs;
+      for (int i = 0; i < fespace->GetNBE(); i++)
+      {
+         // Get the vector degrees of freedom in the boundary element.
+         fespace->GetBdrElementVDofs(i, vdofs);
+         // Set the boundary values to zero.
+         for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
+      }
+      x -= rdm;
    }
-   Array<int> vdofs;
-   for (int i = 0; i < fespace->GetNBE(); i++)
-   {
-      // Get the vector degrees of freedom in the boundary element.
-      fespace->GetBdrElementVDofs(i, vdofs);
-      // Set the boundary values to zero.
-      for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
-   }
-   x -= rdm;
    x.SetTrueVector();
    x.SetFromTrueVector();
 
-   // 9. Save the starting (prior to the optimization) mesh to a file. This
-   //    output can be viewed later using GLVis: "glvis -m perturbed.mesh".
+   // Save the starting (prior to the optimization) mesh to a file. This
+   // output can be viewed later using GLVis: "glvis -m perturbed.mesh".
    {
       ofstream mesh_ofs("perturbed.mesh");
       mesh->Print(mesh_ofs);
    }
 
-   // 10. Store the starting (prior to the optimization) positions.
-   GridFunction x0(fespace);
-   x0 = x;
+   // Store the starting (prior to the optimization) positions.
+   GridFunction x0(x);
 
-   // 11. Form the integrator that uses the chosen metric and target.
+   // Form the integrator that uses the chosen metric and target.
    real_t min_detJ = -0.1;
    TMOP_QualityMetric *metric = NULL;
    switch (metric_id)
@@ -872,12 +874,12 @@ int main(int argc, char *argv[])
    // normalization factors for these terms as well.
    if (normalization) { tmop_integ->EnableNormalization(x0); }
 
-   // 12. Setup the final NonlinearForm (which defines the integral of interest,
-   //     its first and second derivatives). Here we can use a combination of
-   //     metrics, i.e., optimize the sum of two integrals, where both are
-   //     scaled by used-defined space-dependent weights. Note that there are no
-   //     command-line options for the weights and the type of the second
-   //     metric; one should update those in the code.
+   // Setup the final NonlinearForm (which defines the integral of interest,
+   // its first and second derivatives). Here we can use a combination of
+   // metrics, i.e., optimize the sum of two integrals, where both are
+   // scaled by used-defined space-dependent weights. Note that there are no
+   // command-line options for the weights and the type of the second
+   // metric; one should update those in the code.
    NonlinearForm a(fespace);
    if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    ConstantCoefficient *metric_coeff1 = NULL;
@@ -983,11 +985,11 @@ int main(int argc, char *argv[])
       vis_tmop_metric_s(mesh_poly_deg, *metric, *target_c, *mesh, title, 0);
    }
 
-   // 13. Fix all boundary nodes, or fix only a given component depending on the
-   //     boundary attributes of the given mesh. Attributes 1/2/3 correspond to
-   //     fixed x/y/z components of the node. Attribute 4 corresponds to an
-   //     entirely fixed node. Other boundary attributes do not affect the node
-   //     movement boundary conditions.
+   // Fix all boundary nodes, or fix only a given component depending on the
+   // boundary attributes of the given mesh. Attributes 1/2/3 correspond to
+   // fixed x/y/z components of the node. Attribute 4 corresponds to an
+   // entirely fixed node. Other boundary attributes do not affect the node
+   // movement boundary conditions.
    if (move_bnd == false)
    {
       Array<int> ess_bdr(mesh->bdr_attributes.Max());
@@ -1008,7 +1010,7 @@ int main(int argc, char *argv[])
          if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
          if (attr == 4) { n += nd * dim; }
       }
-      Array<int> ess_vdofs(n);
+      Array<int> vdofs, ess_vdofs(n);
       n = 0;
       for (int i = 0; i < mesh->GetNBE(); i++)
       {
@@ -1136,8 +1138,8 @@ int main(int argc, char *argv[])
    }
    hr_solver.Mult();
 
-   // 15. Save the optimized mesh to a file. This output can be viewed later
-   //     using GLVis: "glvis -m optimized.mesh".
+   // Save the optimized mesh to a file. This output can be viewed later
+   // using GLVis: "glvis -m optimized.mesh".
    {
       ofstream mesh_ofs("optimized.mesh");
       mesh_ofs.precision(14);

miniapps/meshing/pmesh-fitting.cpp

@@ -286,11 +286,11 @@ int main (int argc, char *argv[])
       pmesh->PrintAsSerial(mesh_ofs);
    }
 
-   // 11. Store the starting (prior to the optimization) positions.
+   // Store the starting (prior to the optimization) positions.
    ParGridFunction x0(pfespace);
    x0 = x;
 
-   // 12. Form the integrator that uses the chosen metric and target.
+   // Form the integrator that uses the chosen metric and target.
    TMOP_QualityMetric *metric = NULL;
    switch (metric_id)
    {

miniapps/meshing/pmesh-optimizer.cpp

@@ -339,26 +339,19 @@ int main (int argc, char *argv[])
    Device device(devopt);
    if (myid == 0) { device.Print();}
 
-   // 3. Initialize and refine the starting mesh.
+   // Initialize and refine the starting mesh.
    Mesh *mesh = new Mesh(mesh_file, 1, 1, false);
-   for (int lev = 0; lev < rs_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
+   for (int lev = 0; lev < rs_levels; lev++) { mesh->UniformRefinement(); }
    const int dim = mesh->Dimension();
-
    if (hradaptivity) { mesh->EnsureNCMesh(); }
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
    delete mesh;
-   for (int lev = 0; lev < rp_levels; lev++)
-   {
-      pmesh->UniformRefinement();
-   }
+   for (int lev = 0; lev < rp_levels; lev++) { pmesh->UniformRefinement(); }
 
-   // 4. Define a finite element space on the mesh. Here we use vector finite
-   //    elements which are tensor products of quadratic finite elements. The
-   //    number of components in the vector finite element space is specified by
-   //    the last parameter of the FiniteElementSpace constructor.
+   // Define a finite element space on the mesh. Here we use vector finite
+   // elements which are tensor products of quadratic finite elements. The
+   // number of components in the vector finite element space is specified by
+   // the last parameter of the FiniteElementSpace constructor.
    FiniteElementCollection *fec;
    if (mesh_poly_deg <= 0)
    {
@@ -366,25 +359,26 @@ int main (int argc, char *argv[])
       mesh_poly_deg = 2;
    }
    else { fec = new H1_FECollection(mesh_poly_deg, dim); }
-   ParFiniteElementSpace *pfespace = new ParFiniteElementSpace(pmesh, fec, dim,
-                                                               mesh_node_ordering);
+   //else { fec = new L2_FECollection(mesh_poly_deg, dim); }
+   auto pfespace = new ParFiniteElementSpace(pmesh, fec, dim,
+                                             mesh_node_ordering);
 
-   // 5. Make the mesh curved based on the above finite element space. This
-   //    means that we define the mesh elements through a fespace-based
-   //    transformation of the reference element.
+   // Make the mesh curved based on the above finite element space. This
+   // means that we define the mesh elements through a fespace-based
+   // transformation of the reference element.
    pmesh->SetNodalFESpace(pfespace);
 
-   // 7. Get the mesh nodes (vertices and other degrees of freedom in the finite
-   //    element space) as a finite element grid function in fespace. Note that
-   //    changing x automatically changes the shapes of the mesh elements.
+   // Get the mesh nodes (vertices and other degrees of freedom in the finite
+   // element space) as a finite element grid function in fespace. Note that
+   // changing x automatically changes the shapes of the mesh elements.
    ParGridFunction x(pfespace);
    pmesh->SetNodalGridFunction(&x);
 
-   // 8. Define a vector representing the minimal local mesh size in the mesh
-   //    nodes. We index the nodes using the scalar version of the degrees of
-   //    freedom in pfespace. Note: this is partition-dependent.
+   // Define a vector representing the minimal local mesh size in the mesh
+   // nodes. We index the nodes using the scalar version of the degrees of
+   // freedom in pfespace. Note: this is partition-dependent.
    //
-   //    In addition, compute average mesh size and total volume.
+   // In addition, compute average mesh size and total volume.
    Vector h0(pfespace->GetNDofs());
    h0 = infinity();
    real_t vol_loc = 0.0;
@@ -407,39 +401,43 @@ int main (int argc, char *argv[])
    const real_t small_phys_size = pow(vol_glb, 1.0 / dim) / 100.0;
 
    // 9. Add a random perturbation to the nodes in the interior of the domain.
-   //    We define a random grid function of fespace and make sure that it is
-   //    zero on the boundary and its values are locally of the order of h0.
-   //    The latter is based on the DofToVDof() method which maps the scalar to
-   //    the vector degrees of freedom in pfespace.
-   ParGridFunction rdm(pfespace);
-   rdm.Randomize();
-   rdm -= 0.25; // Shift to random values in [-0.5,0.5].
-   rdm *= jitter;
-   rdm.HostReadWrite();
-   // Scale the random values to be of order of the local mesh size.
-   for (int i = 0; i < pfespace->GetNDofs(); i++)
+   // We define a random grid function of fespace and make sure that it is
+   // zero on the boundary and its values are locally of the order of h0.
+   // The latter is based on the DofToVDof() method which maps the scalar to
+   // the vector degrees of freedom in pfespace.
+   if (jitter > 0.0)
    {
-      for (int d = 0; d < dim; d++)
+      ParGridFunction rdm(pfespace);
+      rdm.Randomize();
+      rdm -= 0.25; // Shift to random values in [-0.5,0.5].
+      rdm *= jitter;
+      rdm.HostReadWrite();
+      // Scale the random values to be of order of the local mesh size.
+      for (int i = 0; i < pfespace->GetNDofs(); i++)
       {
-         rdm(pfespace->DofToVDof(i,d)) *= h0(i);
+         for (int d = 0; d < dim; d++)
+         {
+            rdm(pfespace->DofToVDof(i,d)) *= h0(i);
+         }
       }
+      Array<int> vdofs;
+      for (int i = 0; i < pfespace->GetNBE(); i++)
+      {
+         // Get the vector degrees of freedom in the boundary element.
+         pfespace->GetBdrElementVDofs(i, vdofs);
+         // Set the boundary values to zero.
+         for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
+      }
+      x -= rdm;
    }
-   Array<int> vdofs;
-   for (int i = 0; i < pfespace->GetNBE(); i++)
-   {
-      // Get the vector degrees of freedom in the boundary element.
-      pfespace->GetBdrElementVDofs(i, vdofs);
-      // Set the boundary values to zero.
-      for (int j = 0; j < vdofs.Size(); j++) { rdm(vdofs[j]) = 0.0; }
-   }
-   x -= rdm;
+
    // Set the perturbation of all nodes from the true nodes.
    x.SetTrueVector();
    x.SetFromTrueVector();
 
-   // 10. Save the starting (prior to the optimization) mesh to a file. This
-   //     output can be viewed later using GLVis: "glvis -m perturbed -np
-   //     num_mpi_tasks".
+   // Save the starting (prior to the optimization) mesh to a file. This
+   // output can be viewed later using GLVis: "glvis -m perturbed -np
+   // num_mpi_tasks".
    {
       ostringstream mesh_name;
       mesh_name << "perturbed.mesh";
@@ -448,9 +446,8 @@ int main (int argc, char *argv[])
       pmesh->PrintAsOne(mesh_ofs);
    }
 
-   // 11. Store the starting (prior to the optimization) positions.
-   ParGridFunction x0(pfespace);
-   x0 = x;
+   // Store the starting (prior to the optimization) positions.
+   ParGridFunction x0(x);
 
    // 12. Form the integrator that uses the chosen metric and target.
    real_t min_detJ = -0.1;
@@ -911,12 +908,12 @@ int main (int argc, char *argv[])
    // normalization factors for these terms as well.
    if (normalization) { tmop_integ->ParEnableNormalization(x0); }
 
-   // 13. Setup the final NonlinearForm (which defines the integral of interest,
-   //     its first and second derivatives). Here we can use a combination of
-   //     metrics, i.e., optimize the sum of two integrals, where both are
-   //     scaled by used-defined space-dependent weights.  Note that there are
-   //     no command-line options for the weights and the type of the second
-   //     metric; one should update those in the code.
+   // Setup the final NonlinearForm (which defines the integral of interest,
+   // its first and second derivatives). Here we can use a combination of
+   // metrics, i.e., optimize the sum of two integrals, where both are
+   // scaled by used-defined space-dependent weights.  Note that there are
+   // no command-line options for the weights and the type of the second
+   // metric; one should update those in the code.
    ParNonlinearForm a(pfespace);
    if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    ConstantCoefficient *metric_coeff1 = NULL;
@@ -1030,10 +1027,10 @@ int main (int argc, char *argv[])
       vis_tmop_metric_p(mesh_poly_deg, *metric, *target_c, *pmesh, title, 0);
    }
 
-   // 14. Fix all boundary nodes, or fix only a given component depending on the
-   //     boundary attributes of the given mesh.  Attributes 1/2/3 correspond to
-   //     fixed x/y/z components of the node.  Attribute dim+1 corresponds to
-   //     an entirely fixed node.
+   // Fix all boundary nodes, or fix only a given component depending on the
+   // boundary attributes of the given mesh.  Attributes 1/2/3 correspond to
+   // fixed x/y/z components of the node.  Attribute dim+1 corresponds to
+   // an entirely fixed node.
    if (move_bnd == false)
    {
       Array<int> ess_bdr(pmesh->bdr_attributes.Max());
@@ -1054,7 +1051,7 @@ int main (int argc, char *argv[])
          if (attr == 1 || attr == 2 || attr == 3) { n += nd; }
          if (attr == 4) { n += nd * dim; }
       }
-      Array<int> ess_vdofs(n);
+      Array<int> vdofs, ess_vdofs(n);
       n = 0;
       for (int i = 0; i < pmesh->GetNBE(); i++)
       {
@@ -1184,8 +1181,8 @@ int main (int argc, char *argv[])
    }
    hr_solver.Mult();
 
-   // 16. Save the optimized mesh to a file. This output can be viewed later
-   //     using GLVis: "glvis -m optimized -np num_mpi_tasks".
+   // Save the optimized mesh to a file. This output can be viewed later
+   // using GLVis: "glvis -m optimized -np num_mpi_tasks".
    {
       ostringstream mesh_name;
       mesh_name << "optimized.mesh";