Merge pull request #4983 from dnpham23/master

Methods for extracting boundary edge DOFs from a 3D ParMesh and for synchronizing boundary conditions across processor boundaries
Address review feedback: Fix H1 DOF dedup and single-precision reductions
2026-08-20 15:23:36 -07:00 · 2026-08-19 17:16:02 -07:00 · 2026-08-19 16:57:41 -07:00 · 2026-08-19 15:55:05 -07:00 · 2026-08-19 13:30:32 -07:00 · 2026-08-19 12:15:43 -07:00
38 changed files with 2729 additions and 952 deletions
@@ -53,6 +53,16 @@ Discretization improvements
  ComplexHypreParMatrix::GetSystemMatrix, which previously assumed equal
  trial and test spaces.

+- Added FiniteElementSpace::GetBoundaryLoopEdgeDofs to extract the edge DOFs on
+  the perimeter loop of a set of boundary elements, with a ParFiniteElementSpace
+  overload that reconciles the selection across processor boundaries so the
+  result is partition invariant. This is useful for imposing boundary conditions
+  on boundary edge DOFs.
+
+- Added a MaxAbs reduction to GroupCommunicator that selects the signed value of
+  largest magnitude across a group, keeping its sign. Equal-magnitude ties
+  resolve deterministically to the positive value.
+
 Meshing improvements
 --------------------
 - Added support for nonuniform anisotropic mesh refinement on parallel quad/hex
@@ -599,11 +599,6 @@ if (MFEM_USE_ENZYME)
    set(ENZYME_INCLUDE_DIRS ${ENZYME_DIR}/include)
 endif()

-# GLVis
-if (MFEM_USE_GLVIS AND NOT MFEM_FETCH_GLVIS)
-  find_package(GLVis REQUIRED)
-endif()
-
 # MFEM_TIMER_TYPE
 if (NOT DEFINED MFEM_TIMER_TYPE)
  if (APPLE)
@@ -645,7 +640,7 @@ set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
    NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
    ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS CUDSS CALIPER CODIPACK
    BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ
-    ALGOIM ENZYME GLVIS CUDA::cudart)
+    ALGOIM ENZYME CUDA::cudart)

 # Add all created targets and *_FOUND libraries in the variables TPL_TARGETS and
 # TPL_LIBRARIES, respectively.
@@ -825,11 +820,6 @@ endif()

 set(MFEM_CUSTOM_TARGET_PREFIX CACHE STRING "")

-if (MFEM_USE_GLVIS AND MFEM_FETCH_GLVIS)
-    # needs to be after mfem_add_library(mfem) to be able to depend on it
-    find_package(GLVis REQUIRED)
-endif()
-
 #-------------------------------------------------------------------------------
 # Examples, miniapps, benchmarks and testing
 #-------------------------------------------------------------------------------
@@ -930,10 +920,6 @@ target_include_directories(mfem BEFORE
  PUBLIC
    $<INSTALL_INTERFACE:${INSTALL_INCLUDE_DIR}>)

-if (MFEM_USE_GLVIS AND MFEM_FETCH_GLVIS)
-    target_link_libraries(mfem PUBLIC GLVIS)
-endif()
-
 # The 'install' target will not depend on 'all'.
 # set(CMAKE_SKIP_INSTALL_ALL_DEPENDENCY TRUE)

@@ -626,9 +626,6 @@ MFEM_USE_ENZYME = YES/NO
   config/defaults.mk. For more detailed instructions, see the section "Specific
   options for Enzyme" below.

-MFEM_USE_GLVIS = YES/NO
-   Enables using a GLVis stream directly instead of the socket one.
-
 MFEM_BUILD_TAG = (any value)
   An optional tag to characterize the build. Exported to config/config.mk.
   Can be used to identify the MFEM build from other makefiles.
@@ -1095,7 +1092,6 @@ MFEM_USE_BENCHMARK
 MFEM_USE_PARELAG
 MFEM_USE_TRIBOL
 MFEM_USE_ENZYME
-MFEM_USE_GLVIS

 The following options are CMake specific:

@@ -1106,7 +1102,6 @@ MFEM_FETCH_TPLS      - Enable fetching of all supported third-party libraries.
 MFEM_FETCH_GSLIB     - Enable fetching of gslib.
 MFEM_FETCH_HYPRE     - Enable fetching of hypre.
 MFEM_FETCH_METIS     - Enable fetching of metis.
-MFEM_FETCH_GLVIS     - Enable fetching of GLVis.

 External libraries (CMake):
 ---------------------------
@@ -222,7 +222,4 @@
 // Enable Enzyme for AD
 #cmakedefine MFEM_USE_ENZYME

-// Enable GLVis
-#cmakedefine MFEM_USE_GLVIS
-
 #endif // MFEM_CONFIG_HEADER
@@ -1,135 +0,0 @@
-# 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.
-
-# Defines the following variables:
-#   - GLVIS_FOUND
-#   - GLVIS_LIBRARIES
-#   - GLVIS_INCLUDE_DIRS
-
-if (MFEM_FETCH_GLVIS OR MFEM_FETCH_TPLS)
-    message(STATUS "GLVis: Fetch/ExternalProject")
-
-    get_directory_property(COMPILE_OPTS COMPILE_OPTIONS)
-    string(REPLACE ";" " " COMPILE_CXX_FLAGS "${COMPILE_OPTS}")
-
-    # get_directory_property(COMPILE_DEFINITIONS COMPILE_DEFINITIONS)
-    # string(REPLACE "\"" "\\\"" COMPILE_DEFS_QUOTED_STR "${COMPILE_DEFINITIONS}")
-    # string(REPLACE ";" " -D" COMPILE_DEFS_STR "-D${COMPILE_DEFS_QUOTED_STR}")
-    # string(JOIN " " COMPILE_CXX_FLAGS ${COMPILE_OPTS_STR} ${COMPILE_DEFS_STR})
-
-    cmake_host_system_information(RESULT NCPU QUERY NUMBER_OF_LOGICAL_CORES)
-
-    add_library(GLVIS STATIC IMPORTED)
-
-    include(ExternalProject)
-
-    set(FETCH_DIR ${CMAKE_CURRENT_BINARY_DIR}/fetch)
-    set(FETCH_GLVIS "${FETCH_DIR}/glvis")
-
-    ExternalProject_Add(glvis
-        GIT_REPOSITORY https://github.com/GLVis/glvis.git
-        GIT_TAG stream_sessions
-        GIT_SHALLOW TRUE
-        UPDATE_DISCONNECTED TRUE
-        CMAKE_GENERATOR ${CMAKE_GENERATOR}
-        PREFIX ${FETCH_GLVIS}
-        SOURCE_DIR ${FETCH_GLVIS}/src
-        STAMP_DIR ${FETCH_GLVIS}/stamp
-        BINARY_DIR ${FETCH_GLVIS}/build
-        DEPENDS mfem
-        CMAKE_ARGS
-            # -DCMAKE_VERBOSE_MAKEFILE=ON
-            -DMFEM_DIR=${CMAKE_CURRENT_BINARY_DIR}
-            -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
-            -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
-            -DCMAKE_CXX_FLAGS:STRING=${COMPILE_CXX_FLAGS}
-            -DGLVIS_BUILD_LIB_ONLY=ON
-        BUILD_COMMAND ${CMAKE_COMMAND}
-            --build ${FETCH_GLVIS}/build
-            --config $<CONFIG>
-            --target glvis glvis_logo
-            --parallel ${NCPU}
-        BUILD_BYPRODUCTS 
-            ${FETCH_GLVIS}/build/lib/libglvis.a
-            ${FETCH_GLVIS}/build/share/libglvis_logo.a
-        INSTALL_COMMAND "")
-        
-    set_target_properties(GLVIS PROPERTIES
-        IMPORTED_LOCATION ${FETCH_GLVIS}/build/lib/libglvis.a)
-
-    find_package(OpenGL REQUIRED)
-    find_package(GLEW   REQUIRED)
-    find_package(SDL2   REQUIRED)
-    find_package(PNG    REQUIRED)
-    find_package(Freetype REQUIRED)
-    find_package(Fontconfig REQUIRED)
-    if(APPLE)
-        find_library(COCOA_LIBRARY Cocoa)
-    endif()
-    target_link_libraries(GLVIS INTERFACE
-        ${FETCH_GLVIS}/build/share/libglvis_logo.a
-        OpenGL::GL
-        GLEW::GLEW
-        SDL2::SDL2
-        PNG::PNG
-        Freetype::Freetype
-        Fontconfig::Fontconfig)
-    if(APPLE)
-        target_link_libraries(GLVIS INTERFACE ${COCOA_LIBRARY})
-    endif()
-
-    set(GLVIS_FOUND TRUE)
-    return()
-endif()
-
-message(STATUS "[🔵 GLVis 🔵] Find pre-installed package")
-
-include(MfemCmakeUtilities)
-mfem_find_package(GLVis GLVIS GLVIS_DIR
-        "include" "lib/glwindow.hpp"
-        "lib" "build/lib/libglvis.a"
-        "Paths to headers required by GLVis"
-        "Libraries required by GLVis")
-
-if (GLVIS_FOUND)
-        set(GLVIS_INCLUDE_DIRS ${GLVIS_INCLUDE_DIRS}/lib)
-
-        find_library(GLVIS_LOGO_LIBRARY
-                NAMES glvis_logo
-                PATHS ${GLVIS_DIR}/build/share
-                NO_DEFAULT_PATH 
-                REQUIRED)
-        list(APPEND GLVIS_LIBRARIES ${GLVIS_LOGO_LIBRARY})
-
-        find_package(OpenGL REQUIRED)
-        list(APPEND GLVIS_LIBRARIES OpenGL::GL)
-
-        find_package(GLEW REQUIRED)
-        list(APPEND GLVIS_LIBRARIES GLEW::GLEW)
-        
-        find_package(SDL2 REQUIRED)
-        list(APPEND GLVIS_LIBRARIES SDL2::SDL2)
-
-        find_package(PNG REQUIRED)
-        list(APPEND GLVIS_LIBRARIES PNG::PNG)
-
-        find_package(Freetype REQUIRED)
-        list(APPEND GLVIS_LIBRARIES Freetype::Freetype)
-
-        find_package(Fontconfig REQUIRED)
-        list(APPEND GLVIS_LIBRARIES Fontconfig::Fontconfig)
-
-        find_library(COCOA_LIBRARY Cocoa)
-        list(APPEND GLVIS_LIBRARIES ${COCOA_LIBRARY})
-endif()
-        
-message(STATUS "GLVIS_INCLUDE_DIRS: ${GLVIS_INCLUDE_DIRS}")
-message(STATUS "GLVIS_LIBRARIES: ${GLVIS_LIBRARIES}")
@@ -884,7 +884,7 @@ function(mfem_export_mk_files)
      MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
      MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG
      MFEM_USE_TRIBOL MFEM_USE_MOONOLITH MFEM_USE_ALGOIM MFEM_USE_ENZYME
-      MFEM_USE_HDF5 MFEM_USE_GLVIS)
+      MFEM_USE_HDF5)
  foreach(var ${CONFIG_MK_BOOL_VARS})
    if (${var})
      set(${var} YES)
@@ -18,7 +18,7 @@
 #define MFEM_CONFIG_HPP

 #ifdef MFEM_CONFIG_FILE
-#include MFEM_CONFIG_FILE // IWYU pragma: export
+#include MFEM_CONFIG_FILE
 #else
 #include "_config.hpp"
 #endif
@@ -222,7 +222,4 @@
 // Enable the Enzyme LLVM plugin
 // #define MFEM_USE_ENZYME

-// Enable GLVis.
-// #define MFEM_USE_GLVIS
-
 #endif // MFEM_CONFIG_HEADER
@@ -72,7 +72,6 @@ MFEM_USE_BENCHMARK     = @MFEM_USE_BENCHMARK@
 MFEM_USE_PARELAG       = @MFEM_USE_PARELAG@
 MFEM_USE_TRIBOL        = @MFEM_USE_TRIBOL@
 MFEM_USE_ENZYME        = @MFEM_USE_ENZYME@
-MFEM_USE_GLVIS         = @MFEM_USE_GLVIS@

 # Compiler, compile options, and link options
 MFEM_CXX       = @MFEM_CXX@
@@ -72,7 +72,6 @@ option(MFEM_USE_BENCHMARK "Enable Google Benchmark" OFF)
 option(MFEM_USE_PARELAG "Enable ParELAG" OFF)
 option(MFEM_USE_TRIBOL "Enable Tribol" OFF)
 option(MFEM_USE_ENZYME "Enable Enzyme" OFF)
-option(MFEM_USE_GLVIS "Enable GLVis" OFF)

 # Optional overrides for autodetected MPIEXEC and MPIEXEC_NUMPROC_FLAG
 # set(MFEM_MPIEXEC "mpirun" CACHE STRING "Command for running MPI tests")
@@ -97,7 +96,6 @@ option(MFEM_FETCH_TPLS "Enable fetching of all supported third-party libraries"
 option(MFEM_FETCH_GSLIB "Enable fetching of GSLIB" OFF)
 option(MFEM_FETCH_HYPRE "Enable fetching of hypre" OFF)
 option(MFEM_FETCH_METIS "Enable fetching of METIS" OFF)
-option(MFEM_FETCH_GLVIS "Enable fetching of GLVis" OFF)

 # Setting CXX/MPICXX on the command line or in user.cmake will overwrite the
 # autodetected C++ compiler.
@@ -280,8 +278,6 @@ set(Tribol_REQUIRED_PACKAGES "Axom/core/mint/slam/slic" CACHE STRING

 set(ENZYME_DIR "${MFEM_DIR}/../enzyme" CACHE PATH "Path to Enzyme")

-set(GLVIS_DIR "${MFEM_DIR}/../glvis" CACHE PATH "Path to GLVis")
-
 set(BLAS_INCLUDE_DIRS "" CACHE STRING "Path to BLAS headers.")
 set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")
 set(LAPACK_INCLUDE_DIRS "" CACHE STRING "Path to LAPACK headers.")
@@ -207,7 +207,6 @@ MFEM_USE_BENCHMARK     = NO
 MFEM_USE_PARELAG       = NO
 MFEM_USE_TRIBOL        = NO
 MFEM_USE_ENZYME        = NO
-MFEM_USE_GLVIS         = NO

 # Process MFEM_PRECISION -> MFEM_USE_SINGLE, MFEM_USE_DOUBLE
 ifneq ($(filter double Double DOUBLE,$(MFEM_PRECISION)),)
@@ -674,9 +673,6 @@ endif
 ENZYME_OPT = -fplugin=$(ENZYME_PLUGIN)
 ENZYME_LIB =

-# GLVis configuration
-include $(GLVIS_MK)
-
 # If YES, enable some informational messages
 VERBOSE = NO

@@ -1,116 +0,0 @@
-# GLVis library - Adapted from GLVis' makefile
-
-# Macro that searches for a file in a list of directories returning the first
-# directory that contains the file.
-# $(1) - the file to search for
-# $(2) - list of directories to search
-define find_dir
-$(patsubst %/$(1),%,$(firstword $(wildcard $(foreach d,$(2),$(d)/$(1)))))
-endef
-
-# Macro to find the proper library sub-directory, 'lib64' or 'lib', given a
-# tentative prefix and a library name. Returns empty path if prefix is empty,
-# '/usr', or the library is not found.
-# $(1) - the prefix to search, e.g. $(SDL_DIR)
-# $(2) - library name without 'lib' prefix, e.g. 'SDL2'
-define dir2lib
-$(if $(filter-out /usr,$(1)),$(patsubst %/,%,$(dir $(firstword $(wildcard\
- $(1)/lib64/lib$(2).* $(1)/lib/lib$(2).*)))))
-endef
-
-BREW_PREFIX := $(if $(NOTMAC),,$(shell brew --prefix 2> /dev/null))
-
-FREETYPE_SEARCH_PATHS = $(BREW_PREFIX) /usr /opt/X11
-FREETYPE_SEARCH_FILE = include/freetype2/ft2build.h
-FREETYPE_DIR = $(call find_dir,$(FREETYPE_SEARCH_FILE),$(FREETYPE_SEARCH_PATHS))
-FREETYPE_LIB_DIR = $(call dir2lib,$(FREETYPE_DIR),freetype)
-FREETYPE_LIBS = -lfreetype -lfontconfig
-
-# If GLEW is in /usr, there's no need to add search paths
-GLEW_SEARCH_PATHS = /usr/local $(BREW_PREFIX) $(abspath ../glew)
-GLEW_SEARCH_FILE = include/GL/glew.h
-GLEW_DIR ?= $(call find_dir,$(GLEW_SEARCH_FILE),$(GLEW_SEARCH_PATHS))
-GLEW_LIB_DIR = $(call dir2lib,$(GLEW_DIR),GLEW)
-GLEW_LIBS = -lGLEW
-
-# If SDL is in /usr, there's no need to add search paths
-SDL_SEARCH_PATHS := /usr/local $(BREW_PREFIX) $(abspath ../SDL2)
-SDL_SEARCH_FILE = include/SDL2/SDL.h
-SDL_DIR ?= $(call find_dir,$(SDL_SEARCH_FILE),$(SDL_SEARCH_PATHS))
-SDL_LIB_DIR = $(call dir2lib,$(SDL_DIR),SDL2)
-SDL_LIBS = -lSDL2
-
-# If GLM is in /usr/include, there's no need to add search paths
-GLM_SEARCH_PATHS = /usr/local/include \
- $(if $(BREW_PREFIX),$(BREW_PREFIX)/include) $(abspath ../glm)
-GLM_SEARCH_FILE = glm/glm.hpp
-GLM_DIR ?= $(call find_dir,$(GLM_SEARCH_FILE),$(GLM_SEARCH_PATHS))
-
-# If OpenGL is in /usr, there's no need to add search paths
-OPENGL_SEARCH_PATHS = /usr/local /opt/local
-OPENGL_SEARCH_FILE = include/GL/gl.h
-OPENGL_DIR ?= $(call find_dir,$(OPENGL_SEARCH_FILE),$(OPENGL_SEARCH_PATHS))
-OPENGL_LIB_DIR = $(if $(NOTMAC),$(call dir2lib,$(OPENGL_DIR),GL))
-OPENGL_LIBS = $(if $(NOTMAC),-lGL,-framework OpenGL -framework Cocoa)
-
-# Regarding -DGLEW_NO_GLU, see https://github.com/nigels-com/glew/issues/192
-GL_OPTS ?= $(if $(FREETYPE_DIR),-I$(FREETYPE_DIR)/include/freetype2) \
- $(if $(SDL_DIR),-I$(SDL_DIR)/include) \
- $(if $(GLEW_DIR),-I$(GLEW_DIR)/include) -DGLEW_NO_GLU \
- $(if $(GLM_DIR),-I$(GLM_DIR)) \
- $(if $(OPENGL_DIR),-I$(OPENGL_DIR)/include)
-
-rpath=-Wl,-rpath,
-GL_LIBS ?= $(if $(FREETYPE_LIB_DIR),-L$(FREETYPE_LIB_DIR)) \
- $(if $(SDL_LIB_DIR),-L$(SDL_LIB_DIR) $(rpath)$(SDL_LIB_DIR)) \
- $(if $(NOTMAC),$(if $(OPENGL_LIB_DIR),-L$(OPENGL_LIB_DIR) \
-   $(rpath)$(OPENGL_LIB_DIR))) \
- $(if $(GLEW_LIB_DIR),-L$(GLEW_LIB_DIR) $(rpath)$(GLEW_LIB_DIR)) \
- $(FREETYPE_LIBS) $(SDL_LIBS) $(GLEW_LIBS)
-
-GLVIS_FLAGS += $(GL_OPTS)
-GLVIS_LIBS  += $(GL_LIBS)
-
-# Take screenshots internally with libtiff, libpng, or sdl2?
-GLVIS_USE_LIBTIFF ?= NO
-GLVIS_USE_LIBPNG  ?= YES
-TIFF_OPTS = -DGLVIS_USE_LIBTIFF -I/sw/include
-TIFF_LIBS = -L/sw/lib -ltiff
-PNG_OPTS = -DGLVIS_USE_LIBPNG
-PNG_LIBS = -lpng
-ifeq ($(GLVIS_USE_LIBTIFF),YES)
-   GLVIS_FLAGS += $(TIFF_OPTS)
-   GLVIS_LIBS  += $(TIFF_LIBS)
-else ifeq ($(GLVIS_USE_LIBPNG),YES)
-   GLVIS_FLAGS += $(PNG_OPTS)
-   GLVIS_LIBS  += $(PNG_LIBS)
-else
-   # no flag --> SDL screenshots
-endif
-
-# EGL headless rendering
-GLVIS_USE_EGL ?= NO
-EGL_OPTS = -DGLVIS_USE_EGL
-EGL_LIBS = -lEGL
-ifeq ($(GLVIS_USE_EGL),YES)
-	GLVIS_FLAGS += $(EGL_OPTS)
-	GLVIS_LIBS  += $(EGL_LIBS)
-endif
-
-# CGL headless rendering
-GLVIS_USE_CGL ?= $(if $(NOTMAC),NO,YES)
-CGL_OPTS = -DGLVIS_USE_CGL
-ifeq ($(GLVIS_USE_CGL),YES)
-	GLVIS_FLAGS += $(CGL_OPTS)
-endif
-
-PTHREAD_LIB = -lpthread
-
-GLVIS_LIBS += $(PTHREAD_LIB)
-GLVIS_LIBS += $(if $(NOTMAC),-lmfem)
-GLVIS_LIBS := $(sort $(GLVIS_LIBS))
-GLVIS_LIBS += $(OPENGL_LIBS)
-
-GLVIS_DIR = @MFEM_DIR@/../glvis
-GLVIS_OPT = 
-GLVIS_LIB = -L$(GLVIS_DIR)/lib -lglvis $(GLVIS_LIBS)
@@ -0,0 +1,131 @@
+// Define the cube sizes
+L_outer = 1.0;
+L_inner = 0.5;
+
+// Set mesh size and algorithm
+mesh_size = 0.4;
+Mesh.Algorithm3D = 1; // Delaunay algorithm for 3D mesh
+Mesh.CharacteristicLengthFactor = 1.0;
+Mesh.MshFileVersion = 2.2;
+
+// Define center point for concentric cubes
+cx = 0.5;
+cy = 0.5;
+cz = 0.5;
+
+// Define the points (vertices of the outer cube)
+Point(1) = {cx-L_outer/2, cy-L_outer/2, cz-L_outer/2, mesh_size};
+Point(2) = {cx+L_outer/2, cy-L_outer/2, cz-L_outer/2, mesh_size};
+Point(3) = {cx+L_outer/2, cy+L_outer/2, cz-L_outer/2, mesh_size};
+Point(4) = {cx-L_outer/2, cy+L_outer/2, cz-L_outer/2, mesh_size};
+Point(5) = {cx-L_outer/2, cy-L_outer/2, cz+L_outer/2, mesh_size};
+Point(6) = {cx+L_outer/2, cy-L_outer/2, cz+L_outer/2, mesh_size};
+Point(7) = {cx+L_outer/2, cy+L_outer/2, cz+L_outer/2, mesh_size};
+Point(8) = {cx-L_outer/2, cy+L_outer/2, cz+L_outer/2, mesh_size};
+
+// Define the points (vertices of the inner cube)
+Point(9) = {cx-L_inner/2, cy-L_inner/2, cz-L_inner/2, mesh_size};
+Point(10) = {cx+L_inner/2, cy-L_inner/2, cz-L_inner/2, mesh_size};
+Point(11) = {cx+L_inner/2, cy+L_inner/2, cz-L_inner/2, mesh_size};
+Point(12) = {cx-L_inner/2, cy+L_inner/2, cz-L_inner/2, mesh_size};
+Point(13) = {cx-L_inner/2, cy-L_inner/2, cz+L_inner/2, mesh_size};
+Point(14) = {cx+L_inner/2, cy-L_inner/2, cz+L_inner/2, mesh_size};
+Point(15) = {cx+L_inner/2, cy+L_inner/2, cz+L_inner/2, mesh_size};
+Point(16) = {cx-L_inner/2, cy+L_inner/2, cz+L_inner/2, mesh_size};
+
+// Define the lines (edges of the outer cube)
+Line(1) = {1, 2};
+Line(2) = {2, 3};
+Line(3) = {3, 4};
+Line(4) = {4, 1};
+Line(5) = {5, 6};
+Line(6) = {6, 7};
+Line(7) = {7, 8};
+Line(8) = {8, 5};
+Line(9) = {1, 5};
+Line(10) = {2, 6};
+Line(11) = {3, 7};
+Line(12) = {4, 8};
+
+// Define the lines (edges of the inner cube)
+Line(13) = {9, 10};
+Line(14) = {10, 11};
+Line(15) = {11, 12};
+Line(16) = {12, 9};
+Line(17) = {13, 14};
+Line(18) = {14, 15};
+Line(19) = {15, 16};
+Line(20) = {16, 13};
+Line(21) = {9, 13};
+Line(22) = {10, 14};
+Line(23) = {11, 15};
+Line(24) = {12, 16};
+
+// Define the surfaces (faces of the outer cube)
+Line Loop(1) = {1, 2, 3, 4};
+Plane Surface(1) = {1};
+
+Line Loop(2) = {5, 6, 7, 8};
+Plane Surface(2) = {2};
+
+Line Loop(3) = {9, 5, -10, -1};
+Plane Surface(3) = {3};
+
+Line Loop(4) = {10, 6, -11, -2};
+Plane Surface(4) = {4};
+
+Line Loop(5) = {11, 7, -12, -3};
+Plane Surface(5) = {5};
+
+Line Loop(6) = {12, 8, -9, -4};
+Plane Surface(6) = {6};
+
+// Define the surfaces (faces of the inner cube)
+Line Loop(7) = {13, 14, 15, 16};
+Plane Surface(7) = {7};
+
+Line Loop(8) = {17, 18, 19, 20};
+Plane Surface(8) = {8};
+
+Line Loop(9) = {21, 17, -22, -13};
+Plane Surface(9) = {9};
+
+Line Loop(10) = {22, 18, -23, -14};
+Plane Surface(10) = {10};
+
+Line Loop(11) = {23, 19, -24, -15};
+Plane Surface(11) = {11};
+
+Line Loop(12) = {24, 20, -21, -16};
+Plane Surface(12) = {12};
+
+// Define the volumes
+Surface Loop(1) = {1, 2, 3, 4, 5, 6};
+Surface Loop(2) = {7, 8, 9, 10, 11, 12};
+Volume(1) = {1, 2}; // Outer volume with inner hole
+Volume(2) = {2};    // Inner volume
+
+// Assign physical groups
+Physical Volume(1) = {1}; // Outer volume
+Physical Volume(2) = {2}; // Inner volume
+
+// Outer cube surfaces
+Physical Surface(1) = {1}; // Outer bottom
+Physical Surface(2) = {2}; // Outer top
+Physical Surface(3) = {3}; // Outer front
+Physical Surface(4) = {4}; // Outer right
+Physical Surface(5) = {5}; // Outer back
+Physical Surface(6) = {6}; // Outer left
+
+// Inner cube surfaces
+Physical Surface(7) = {7};  // Inner bottom (-xy)
+Physical Surface(8) = {8};  // Inner top (+xy)
+Physical Surface(9) = {9};  // Inner front (-xz)
+Physical Surface(10) = {10}; // Inner right (+yz)
+Physical Surface(11) = {11}; // Inner back (+xz)
+Physical Surface(12) = {12}; // Inner left (-yz)
+
+// Mesh control
+Mesh.OptimizeNetgen = 1;
+Mesh.Optimize = 1;
+Mesh.ElementOrder = 1;
@@ -0,0 +1,907 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$Nodes
+138
+1 0 0 0
+2 1 0 0
+3 1 1 0
+4 0 1 0
+5 0 0 1
+6 1 0 1
+7 1 1 1
+8 0 1 1
+9 0.25 0.25 0.25
+10 0.75 0.25 0.25
+11 0.75 0.75 0.25
+12 0.25 0.75 0.25
+13 0.25 0.25 0.75
+14 0.75 0.25 0.75
+15 0.75 0.75 0.75
+16 0.25 0.75 0.75
+17 0.3333333333325025 0 0
+18 0.6666666666657889 0 0
+19 1 0.3333333333325025 0
+20 1 0.6666666666657889 0
+21 0.6666666666675911 1 0
+22 0.3333333333347203 1 0
+23 0 0.6666666666675911 0
+24 0 0.3333333333347203 0
+25 0.3333333333325025 0 1
+26 0.6666666666657889 0 1
+27 1 0.3333333333325025 1
+28 1 0.6666666666657889 1
+29 0.6666666666675911 1 1
+30 0.3333333333347203 1 1
+31 0 0.6666666666675911 1
+32 0 0.3333333333347203 1
+33 0 0 0.3333333333325025
+34 0 0 0.6666666666657889
+35 1 0 0.3333333333325025
+36 1 0 0.6666666666657889
+37 1 1 0.3333333333325025
+38 1 1 0.6666666666657889
+39 0 1 0.3333333333325025
+40 0 1 0.6666666666657889
+41 0.5000000000003468 0.25 0.25
+42 0.75 0.5000000000003468 0.25
+43 0.5000000000013763 0.75 0.25
+44 0.25 0.5000000000013763 0.25
+45 0.5000000000003468 0.25 0.75
+46 0.75 0.5000000000003468 0.75
+47 0.5000000000013763 0.75 0.75
+48 0.25 0.5000000000013763 0.75
+49 0.25 0.25 0.5000000000003468
+50 0.75 0.25 0.5000000000003468
+51 0.75 0.75 0.5000000000003468
+52 0.25 0.75 0.5000000000003468
+53 0.7113248654055673 0.4999999999991457 0
+54 0.2886751345942123 0.5000000000011557 0
+55 0.5000000000006117 0.7525600817161773 0
+56 0.4999999999993867 0.2474399182839603 0
+57 0.2423197548524782 0.7576802451481532 0
+58 0.757680245147464 0.2423197548520695 0
+59 0.2423197548507857 0.2423197548513912 0
+60 0.7576802451491019 0.7576802451486099 0
+61 0.7113248654055673 0.4999999999991457 1
+62 0.2886751345942123 0.5000000000011557 1
+63 0.5000000000006117 0.7525600817161773 1
+64 0.4999999999993867 0.2474399182839603 1
+65 0.2423197548524782 0.7576802451481532 1
+66 0.757680245147464 0.2423197548520695 1
+67 0.2423197548507857 0.2423197548513912 1
+68 0.7576802451491019 0.7576802451486099 1
+69 0.4999999999993203 0 0.301447615129799
+70 0.4999999999992795 0 0.7028666213189801
+71 0.7525600817158393 0 0.5007190394076877
+72 0.2474399182836191 0 0.5007190394076877
+73 0.7576802451479793 0 0.7576802451479793
+74 0.2423197548517962 0 0.7576802451477375
+75 0.7576802451484569 0 0.2423197548510767
+76 0.2423197548513188 0 0.2423197548513187
+77 1 0.4999999999993203 0.301447615129799
+78 1 0.4999999999992795 0.7028666213189801
+79 1 0.7525600817158394 0.5007190394076877
+80 1 0.2474399182836191 0.5007190394076877
+81 1 0.7576802451479794 0.7576802451479794
+82 1 0.2423197548517962 0.7576802451477376
+83 1 0.7576802451484569 0.2423197548510768
+84 1 0.2423197548513188 0.2423197548513188
+85 0.5000000000008327 1 0.3014476151298047
+86 0.500000000000961 1 0.7028666213191928
+87 0.2474399182842484 1 0.5007190394077241
+88 0.7525600817164384 1 0.5007190394078933
+89 0.2423197548520873 1 0.7576802451480517
+90 0.7576802451481496 1 0.2423197548518761
+91 0.2423197548516099 1 0.2423197548510044
+92 0.7576802451486874 1 0.7576802451481952
+93 0 0.5000000000008327 0.3014476151298047
+94 0 0.500000000000961 0.7028666213191928
+95 0 0.2474399182842484 0.5007190394077241
+96 0 0.7525600817164384 0.5007190394078933
+97 0 0.2423197548520873 0.7576802451480517
+98 0 0.7576802451481496 0.2423197548518761
+99 0 0.2423197548516099 0.2423197548510044
+100 0 0.7576802451486874 0.7576802451481952
+101 0.3968750000003409 0.603125000000244 0.25
+102 0.4374999999998713 0.4375000000001287 0.25
+103 0.5739583333335919 0.5718750000001767 0.25
+104 0.6093749999999631 0.3906250000003402 0.25
+105 0.3968750000003409 0.603125000000244 0.75
+106 0.4374999999998713 0.4375000000001287 0.75
+107 0.5739583333335919 0.5718750000001767 0.75
+108 0.6093749999999631 0.3906250000003402 0.75
+109 0.3806942419826734 0.25 0.3806942419826734
+110 0.5625000000001735 0.25 0.4375000000000001
+111 0.4254282069971791 0.25 0.5712615403304835
+112 0.6093749999998808 0.25 0.6093749999998808
+113 0.75 0.3806942419826734 0.3806942419826734
+114 0.75 0.5625000000001735 0.4375000000000001
+115 0.75 0.4254282069971791 0.5712615403304835
+116 0.75 0.6093749999998808 0.6093749999998808
+117 0.3968750000004991 0.75 0.3968750000002804
+118 0.4375000000000001 0.75 0.5625000000004308
+119 0.5739583333336153 0.75 0.4281250000000707
+120 0.6093749999998166 0.75 0.6093749999993661
+121 0.25 0.3968750000004991 0.3968750000002804
+122 0.25 0.4375000000000001 0.5625000000004308
+123 0.25 0.5739583333336153 0.4281250000000707
+124 0.25 0.6093749999998166 0.6093749999993661
+125 0.4962939304035875 0.5214350017087855 0.4925553109323813
+126 0.3432581985549767 0.6471275530923523 0.3554206606168006
+127 0.6442168181713744 0.5929232373214773 0.3593785632300704
+128 0.625174517421737 0.3491579444372839 0.3604636129462503
+129 0.6130544111091688 0.6576364353993367 0.5046720648819553
+130 0.4281518698369243 0.3632662294430484 0.3548726263205102
+131 0.3639383531355198 0.3520949250221278 0.4978389662613994
+132 0.629585530087249 0.3489878162230438 0.5124654846339122
+133 0.3710853378652663 0.6517586292698121 0.6382643241302075
+134 0.5917018263727056 0.6522525456211955 0.6390764961119443
+135 0.3530810314228338 0.4582477062424107 0.6430239639637545
+136 0.6571010904289998 0.5300774811423468 0.6603970500567977
+137 0.6484596018596915 0.361399676127967 0.6360267588157177
+138 0.4782020887035478 0.3534611476388013 0.6141275027013793
+$EndNodes
+$Elements
+760
+1 2 2 1 1 1 17 59
+2 2 2 1 1 24 1 59
+3 2 2 1 1 18 2 58
+4 2 2 1 1 2 19 58
+5 2 2 1 1 20 3 60
+6 2 2 1 1 3 21 60
+7 2 2 1 1 22 4 57
+8 2 2 1 1 4 23 57
+9 2 2 1 1 17 18 56
+10 2 2 1 1 17 56 59
+11 2 2 1 1 56 18 58
+12 2 2 1 1 19 20 53
+13 2 2 1 1 19 53 58
+14 2 2 1 1 53 20 60
+15 2 2 1 1 21 22 55
+16 2 2 1 1 21 55 60
+17 2 2 1 1 55 22 57
+18 2 2 1 1 23 24 54
+19 2 2 1 1 23 54 57
+20 2 2 1 1 54 24 59
+21 2 2 1 1 54 53 55
+22 2 2 1 1 53 54 56
+23 2 2 1 1 55 53 60
+24 2 2 1 1 53 56 58
+25 2 2 1 1 54 55 57
+26 2 2 1 1 56 54 59
+27 2 2 2 2 5 25 67
+28 2 2 2 2 32 5 67
+29 2 2 2 2 26 6 66
+30 2 2 2 2 6 27 66
+31 2 2 2 2 28 7 68
+32 2 2 2 2 7 29 68
+33 2 2 2 2 30 8 65
+34 2 2 2 2 8 31 65
+35 2 2 2 2 25 26 64
+36 2 2 2 2 25 64 67
+37 2 2 2 2 64 26 66
+38 2 2 2 2 27 28 61
+39 2 2 2 2 27 61 66
+40 2 2 2 2 61 28 68
+41 2 2 2 2 29 30 63
+42 2 2 2 2 29 63 68
+43 2 2 2 2 63 30 65
+44 2 2 2 2 31 32 62
+45 2 2 2 2 31 62 65
+46 2 2 2 2 62 32 67
+47 2 2 2 2 62 61 63
+48 2 2 2 2 61 62 64
+49 2 2 2 2 63 61 68
+50 2 2 2 2 61 64 66
+51 2 2 2 2 62 63 65
+52 2 2 2 2 64 62 67
+53 2 2 3 3 17 1 76
+54 2 2 3 3 1 33 76
+55 2 2 3 3 2 18 75
+56 2 2 3 3 35 2 75
+57 2 2 3 3 5 25 74
+58 2 2 3 3 34 5 74
+59 2 2 3 3 26 6 73
+60 2 2 3 3 6 36 73
+61 2 2 3 3 18 17 69
+62 2 2 3 3 69 17 76
+63 2 2 3 3 18 69 75
+64 2 2 3 3 25 26 70
+65 2 2 3 3 25 70 74
+66 2 2 3 3 70 26 73
+67 2 2 3 3 33 34 72
+68 2 2 3 3 33 72 76
+69 2 2 3 3 72 34 74
+70 2 2 3 3 36 35 71
+71 2 2 3 3 71 35 75
+72 2 2 3 3 36 71 73
+73 2 2 3 3 69 70 71
+74 2 2 3 3 70 69 72
+75 2 2 3 3 69 71 75
+76 2 2 3 3 72 69 76
+77 2 2 3 3 71 70 73
+78 2 2 3 3 70 72 74
+79 2 2 4 4 19 2 84
+80 2 2 4 4 2 35 84
+81 2 2 4 4 3 20 83
+82 2 2 4 4 37 3 83
+83 2 2 4 4 6 27 82
+84 2 2 4 4 36 6 82
+85 2 2 4 4 28 7 81
+86 2 2 4 4 7 38 81
+87 2 2 4 4 20 19 77
+88 2 2 4 4 77 19 84
+89 2 2 4 4 20 77 83
+90 2 2 4 4 27 28 78
+91 2 2 4 4 27 78 82
+92 2 2 4 4 78 28 81
+93 2 2 4 4 35 36 80
+94 2 2 4 4 35 80 84
+95 2 2 4 4 80 36 82
+96 2 2 4 4 38 37 79
+97 2 2 4 4 79 37 83
+98 2 2 4 4 38 79 81
+99 2 2 4 4 77 78 79
+100 2 2 4 4 78 77 80
+101 2 2 4 4 77 79 83
+102 2 2 4 4 80 77 84
+103 2 2 4 4 79 78 81
+104 2 2 4 4 78 80 82
+105 2 2 5 5 21 3 90
+106 2 2 5 5 3 37 90
+107 2 2 5 5 4 22 91
+108 2 2 5 5 39 4 91
+109 2 2 5 5 7 29 92
+110 2 2 5 5 38 7 92
+111 2 2 5 5 30 8 89
+112 2 2 5 5 8 40 89
+113 2 2 5 5 22 21 85
+114 2 2 5 5 85 21 90
+115 2 2 5 5 22 85 91
+116 2 2 5 5 29 30 86
+117 2 2 5 5 29 86 92
+118 2 2 5 5 86 30 89
+119 2 2 5 5 37 38 88
+120 2 2 5 5 37 88 90
+121 2 2 5 5 88 38 92
+122 2 2 5 5 40 39 87
+123 2 2 5 5 87 39 91
+124 2 2 5 5 40 87 89
+125 2 2 5 5 85 86 87
+126 2 2 5 5 86 85 88
+127 2 2 5 5 85 87 91
+128 2 2 5 5 88 85 90
+129 2 2 5 5 87 86 89
+130 2 2 5 5 86 88 92
+131 2 2 6 6 1 24 99
+132 2 2 6 6 33 1 99
+133 2 2 6 6 23 4 98
+134 2 2 6 6 4 39 98
+135 2 2 6 6 32 5 97
+136 2 2 6 6 5 34 97
+137 2 2 6 6 8 31 100
+138 2 2 6 6 40 8 100
+139 2 2 6 6 24 23 93
+140 2 2 6 6 93 23 98
+141 2 2 6 6 24 93 99
+142 2 2 6 6 31 32 94
+143 2 2 6 6 31 94 100
+144 2 2 6 6 94 32 97
+145 2 2 6 6 34 33 95
+146 2 2 6 6 95 33 99
+147 2 2 6 6 34 95 97
+148 2 2 6 6 39 40 96
+149 2 2 6 6 39 96 98
+150 2 2 6 6 96 40 100
+151 2 2 6 6 93 94 95
+152 2 2 6 6 94 93 96
+153 2 2 6 6 93 95 99
+154 2 2 6 6 96 93 98
+155 2 2 6 6 95 94 97
+156 2 2 6 6 94 96 100
+157 2 2 7 7 9 41 102
+158 2 2 7 7 44 9 102
+159 2 2 7 7 41 10 104
+160 2 2 7 7 10 42 104
+161 2 2 7 7 42 11 103
+162 2 2 7 7 11 43 103
+163 2 2 7 7 43 12 101
+164 2 2 7 7 12 44 101
+165 2 2 7 7 102 41 104
+166 2 2 7 7 42 103 104
+167 2 2 7 7 43 101 103
+168 2 2 7 7 101 44 102
+169 2 2 7 7 101 102 103
+170 2 2 7 7 103 102 104
+171 2 2 8 8 13 45 106
+172 2 2 8 8 48 13 106
+173 2 2 8 8 45 14 108
+174 2 2 8 8 14 46 108
+175 2 2 8 8 46 15 107
+176 2 2 8 8 15 47 107
+177 2 2 8 8 47 16 105
+178 2 2 8 8 16 48 105
+179 2 2 8 8 106 45 108
+180 2 2 8 8 46 107 108
+181 2 2 8 8 47 105 107
+182 2 2 8 8 105 48 106
+183 2 2 8 8 105 106 107
+184 2 2 8 8 107 106 108
+185 2 2 9 9 41 9 109
+186 2 2 9 9 9 49 109
+187 2 2 9 9 10 41 110
+188 2 2 9 9 50 10 110
+189 2 2 9 9 13 45 111
+190 2 2 9 9 49 13 111
+191 2 2 9 9 45 14 112
+192 2 2 9 9 14 50 112
+193 2 2 9 9 41 109 110
+194 2 2 9 9 111 45 112
+195 2 2 9 9 109 49 111
+196 2 2 9 9 50 110 112
+197 2 2 9 9 110 109 111
+198 2 2 9 9 110 111 112
+199 2 2 10 10 42 10 113
+200 2 2 10 10 10 50 113
+201 2 2 10 10 11 42 114
+202 2 2 10 10 51 11 114
+203 2 2 10 10 14 46 115
+204 2 2 10 10 50 14 115
+205 2 2 10 10 46 15 116
+206 2 2 10 10 15 51 116
+207 2 2 10 10 42 113 114
+208 2 2 10 10 115 46 116
+209 2 2 10 10 113 50 115
+210 2 2 10 10 51 114 116
+211 2 2 10 10 114 113 115
+212 2 2 10 10 114 115 116
+213 2 2 11 11 43 11 119
+214 2 2 11 11 11 51 119
+215 2 2 11 11 12 43 117
+216 2 2 11 11 52 12 117
+217 2 2 11 11 15 47 120
+218 2 2 11 11 51 15 120
+219 2 2 11 11 47 16 118
+220 2 2 11 11 16 52 118
+221 2 2 11 11 117 43 119
+222 2 2 11 11 47 118 120
+223 2 2 11 11 119 51 120
+224 2 2 11 11 52 117 118
+225 2 2 11 11 118 117 119
+226 2 2 11 11 118 119 120
+227 2 2 12 12 9 44 121
+228 2 2 12 12 49 9 121
+229 2 2 12 12 44 12 123
+230 2 2 12 12 12 52 123
+231 2 2 12 12 48 13 122
+232 2 2 12 12 13 49 122
+233 2 2 12 12 16 48 124
+234 2 2 12 12 52 16 124
+235 2 2 12 12 121 44 123
+236 2 2 12 12 48 122 124
+237 2 2 12 12 49 121 122
+238 2 2 12 12 123 52 124
+239 2 2 12 12 122 121 123
+240 2 2 12 12 122 123 124
+241 4 2 1 1 105 62 106 107
+242 4 2 1 1 102 54 101 103
+243 4 2 1 1 118 86 120 119
+244 4 2 1 1 124 94 123 122
+245 4 2 1 1 52 39 12 96
+246 4 2 1 1 52 12 39 87
+247 4 2 1 1 88 38 15 51
+248 4 2 1 1 79 15 38 51
+249 4 2 1 1 80 14 50 36
+250 4 2 1 1 71 50 14 36
+251 4 2 1 1 120 88 15 51
+252 4 2 1 1 116 15 79 51
+253 4 2 1 1 50 14 112 71
+254 4 2 1 1 52 117 12 87
+255 4 2 1 1 111 69 109 110
+256 4 2 1 1 114 77 115 113
+257 4 2 1 1 124 123 94 96
+258 4 2 1 1 120 86 88 119
+259 4 2 1 1 14 45 64 108
+260 4 2 1 1 103 54 101 55
+261 4 2 1 1 62 105 63 107
+262 4 2 1 1 63 29 15 47
+263 4 2 1 1 14 45 26 64
+264 4 2 1 1 51 11 88 37
+265 4 2 1 1 80 50 10 35
+266 4 2 1 1 49 33 9 72
+267 4 2 1 1 71 10 50 35
+268 4 2 1 1 79 11 51 37
+269 4 2 1 1 49 9 33 95
+270 4 2 1 1 12 123 52 96
+271 4 2 1 1 96 16 52 40
+272 4 2 1 1 49 13 34 72
+273 4 2 1 1 87 52 16 40
+274 4 2 1 1 13 49 34 95
+275 4 2 1 1 43 101 12 55
+276 4 2 1 1 43 12 22 55
+277 4 2 1 1 106 61 108 107
+278 4 2 1 1 102 103 104 53
+279 4 2 1 1 50 115 14 80
+280 4 2 1 1 80 10 50 113
+281 4 2 1 1 9 109 49 72
+282 4 2 1 1 47 15 63 107
+283 4 2 1 1 124 16 52 96
+284 4 2 1 1 49 121 9 95
+285 4 2 1 1 15 28 68 46
+286 4 2 1 1 81 28 15 46
+287 4 2 1 1 92 15 29 47
+288 4 2 1 1 27 82 14 46
+289 4 2 1 1 14 66 27 46
+290 4 2 1 1 26 45 14 73
+291 4 2 1 1 78 115 116 114
+292 4 2 1 1 112 70 111 110
+293 4 2 1 1 79 11 114 51
+294 4 2 1 1 71 50 10 110
+295 4 2 1 1 10 41 104 56
+296 4 2 1 1 77 80 115 113
+297 4 2 1 1 69 109 72 111
+298 4 2 1 1 63 16 30 47
+299 4 2 1 1 25 45 13 64
+300 4 2 1 1 52 16 118 87
+301 4 2 1 1 95 13 49 122
+302 4 2 1 1 57 12 23 44
+303 4 2 1 1 119 11 88 51
+304 4 2 1 1 41 17 9 56
+305 4 2 1 1 18 41 10 56
+306 4 2 1 1 21 11 43 55
+307 4 2 1 1 63 105 16 47
+308 4 2 1 1 13 45 106 64
+309 4 2 1 1 9 102 41 56
+310 4 2 1 1 72 49 13 111
+311 4 2 1 1 11 103 43 55
+312 4 2 1 1 70 112 71 110
+313 4 2 1 1 116 79 78 114
+314 4 2 1 1 32 13 48 67
+315 4 2 1 1 16 89 30 47
+316 4 2 1 1 32 48 13 97
+317 4 2 1 1 48 16 100 31
+318 4 2 1 1 121 123 93 122
+319 4 2 1 1 119 117 118 85
+320 4 2 1 1 10 19 42 84
+321 4 2 1 1 17 9 76 41
+322 4 2 1 1 24 9 59 44
+323 4 2 1 1 42 19 10 58
+324 4 2 1 1 10 41 18 75
+325 4 2 1 1 20 83 11 42
+326 4 2 1 1 90 11 43 21
+327 4 2 1 1 53 104 102 56
+328 4 2 1 1 106 61 64 108
+329 4 2 1 1 87 118 117 85
+330 4 2 1 1 95 121 93 122
+331 4 2 1 1 80 50 115 113
+332 4 2 1 1 109 49 72 111
+333 4 2 1 1 78 27 28 46
+334 4 2 1 1 31 94 48 32
+335 4 2 1 1 19 77 20 42
+336 4 2 1 1 23 24 44 93
+337 4 2 1 1 123 124 52 96
+338 4 2 1 1 88 120 119 51
+339 4 2 1 1 91 43 85 117
+340 4 2 1 1 39 98 12 96
+341 4 2 1 1 12 91 39 87
+342 4 2 1 1 15 92 38 88
+343 4 2 1 1 38 81 15 79
+344 4 2 1 1 36 80 14 82
+345 4 2 1 1 14 71 36 73
+346 4 2 1 1 103 102 54 53
+347 4 2 1 1 61 106 62 107
+348 4 2 1 1 86 29 30 47
+349 4 2 1 1 26 45 70 25
+350 4 2 1 1 21 85 22 43
+351 4 2 1 1 94 93 123 122
+352 4 2 1 1 119 118 86 85
+353 4 2 1 1 12 96 93 123
+354 4 2 1 1 115 78 77 114
+355 4 2 1 1 69 111 70 110
+356 4 2 1 1 28 27 61 46
+357 4 2 1 1 48 62 31 32
+358 4 2 1 1 88 15 92 120
+359 4 2 1 1 79 81 15 116
+360 4 2 1 1 71 14 112 73
+361 4 2 1 1 114 116 79 51
+362 4 2 1 1 112 50 71 110
+363 4 2 1 1 91 12 117 87
+364 4 2 1 1 80 115 14 82
+365 4 2 1 1 63 15 29 68
+366 4 2 1 1 64 26 14 66
+367 4 2 1 1 47 63 105 107
+368 4 2 1 1 43 103 101 55
+369 4 2 1 1 35 10 80 84
+370 4 2 1 1 33 76 9 72
+371 4 2 1 1 37 11 88 90
+372 4 2 1 1 35 71 10 75
+373 4 2 1 1 37 79 11 83
+374 4 2 1 1 9 99 33 95
+375 4 2 1 1 68 63 15 107
+376 4 2 1 1 99 44 93 121
+377 4 2 1 1 18 69 41 17
+378 4 2 1 1 10 35 2 84
+379 4 2 1 1 15 7 28 81
+380 4 2 1 1 3 37 90 11
+381 4 2 1 1 108 64 14 66
+382 4 2 1 1 51 38 79 37
+383 4 2 1 1 12 4 23 98
+384 4 2 1 1 12 39 91 4
+385 4 2 1 1 49 33 72 34
+386 4 2 1 1 10 71 110 75
+387 4 2 1 1 114 11 79 83
+388 4 2 1 1 40 96 16 100
+389 4 2 1 1 34 13 74 72
+390 4 2 1 1 16 87 40 89
+391 4 2 1 1 34 97 13 95
+392 4 2 1 1 20 19 42 53
+393 4 2 1 1 54 23 24 44
+394 4 2 1 1 87 16 118 89
+395 4 2 1 1 95 97 13 122
+396 4 2 1 1 57 22 12 55
+397 4 2 1 1 41 102 104 56
+398 4 2 1 1 64 45 106 108
+399 4 2 1 1 101 57 12 55
+400 4 2 1 1 74 13 5 25
+401 4 2 1 1 13 34 74 5
+402 4 2 1 1 8 40 16 100
+403 4 2 1 1 16 31 65 8
+404 4 2 1 1 82 36 6 14
+405 4 2 1 1 102 9 59 56
+406 4 2 1 1 106 67 13 64
+407 4 2 1 1 119 88 11 90
+408 4 2 1 1 95 49 121 122
+409 4 2 1 1 52 118 117 87
+410 4 2 1 1 124 16 96 100
+411 4 2 1 1 80 10 113 84
+412 4 2 1 1 109 9 76 72
+413 4 2 1 1 72 13 70 111
+414 4 2 1 1 99 9 121 95
+415 4 2 1 1 30 65 16 63
+416 4 2 1 1 13 67 25 64
+417 4 2 1 1 11 55 53 103
+418 4 2 1 1 123 12 44 93
+419 4 2 1 1 56 104 10 58
+420 4 2 1 1 10 18 56 58
+421 4 2 1 1 21 60 11 55
+422 4 2 1 1 17 59 9 56
+423 4 2 1 1 16 63 62 105
+424 4 2 1 1 48 105 16 62
+425 4 2 1 1 111 45 13 70
+426 4 2 1 1 119 86 88 85
+427 4 2 1 1 123 93 94 96
+428 4 2 1 1 115 78 80 77
+429 4 2 1 1 111 72 69 70
+430 4 2 1 1 55 103 54 53
+431 4 2 1 1 63 61 62 107
+432 4 2 1 1 14 27 66 6
+433 4 2 1 1 99 9 1 24
+434 4 2 1 1 49 33 34 95
+435 4 2 1 1 96 52 39 40
+436 4 2 1 1 22 12 91 4
+437 4 2 1 1 20 11 83 3
+438 4 2 1 1 29 15 92 7
+439 4 2 1 1 18 10 75 2
+440 4 2 1 1 36 50 80 35
+441 4 2 1 1 18 17 41 56
+442 4 2 1 1 51 88 38 37
+443 4 2 1 1 62 31 16 48
+444 4 2 1 1 62 16 31 65
+445 4 2 1 1 26 45 25 64
+446 4 2 1 1 104 58 42 10
+447 4 2 1 1 29 63 30 47
+448 4 2 1 1 77 78 79 114
+449 4 2 1 1 71 69 70 110
+450 4 2 1 1 42 103 11 53
+451 4 2 1 1 93 94 95 122
+452 4 2 1 1 118 87 86 85
+453 4 2 1 1 61 62 106 64
+454 4 2 1 1 53 102 54 56
+455 4 2 1 1 78 82 27 46
+456 4 2 1 1 94 48 32 97
+457 4 2 1 1 86 30 89 47
+458 4 2 1 1 78 28 81 46
+459 4 2 1 1 26 70 45 73
+460 4 2 1 1 100 94 48 31
+461 4 2 1 1 86 92 29 47
+462 4 2 1 1 76 69 17 41
+463 4 2 1 1 42 19 77 84
+464 4 2 1 1 77 83 20 42
+465 4 2 1 1 18 41 69 75
+466 4 2 1 1 43 85 90 21
+467 4 2 1 1 91 43 117 12
+468 4 2 1 1 48 62 32 67
+469 4 2 1 1 27 66 61 46
+470 4 2 1 1 68 28 61 46
+471 4 2 1 1 54 24 59 44
+472 4 2 1 1 19 42 53 58
+473 4 2 1 1 23 54 57 44
+474 4 2 1 1 12 93 96 98
+475 4 2 1 1 99 44 121 9
+476 4 2 1 1 42 58 104 53
+477 4 2 1 1 21 43 22 55
+478 4 2 1 1 53 20 11 42
+479 4 2 1 1 53 11 20 60
+480 4 2 1 1 106 67 48 13
+481 4 2 1 1 36 71 50 35
+482 4 2 1 1 52 87 39 40
+483 4 2 1 1 48 67 106 62
+484 4 2 1 1 72 70 13 74
+485 4 2 1 1 9 33 99 1
+486 4 2 1 1 3 37 11 83
+487 4 2 1 1 13 34 5 97
+488 4 2 1 1 12 39 4 98
+489 4 2 1 1 14 36 6 73
+490 4 2 1 1 15 38 92 7
+491 4 2 1 1 8 40 89 16
+492 4 2 1 1 75 35 2 10
+493 4 2 1 1 30 16 65 8
+494 4 2 1 1 14 27 6 82
+495 4 2 1 1 26 14 66 6
+496 4 2 1 1 16 31 8 100
+497 4 2 1 1 32 13 67 5
+498 4 2 1 1 13 67 5 25
+499 4 2 1 1 68 7 28 15
+500 4 2 1 1 29 15 7 68
+501 4 2 1 1 57 4 23 12
+502 4 2 1 1 22 12 4 57
+503 4 2 1 1 10 19 84 2
+504 4 2 1 1 9 59 1 24
+505 4 2 1 1 17 9 59 1
+506 4 2 1 1 18 10 2 58
+507 4 2 1 1 11 21 90 3
+508 4 2 1 1 20 11 3 60
+509 4 2 1 1 16 62 63 65
+510 4 2 1 1 11 53 55 60
+511 4 2 1 1 78 81 79 116
+512 4 2 1 1 71 112 70 73
+513 4 2 1 1 86 88 92 120
+514 4 2 1 1 96 94 124 100
+515 4 2 1 1 94 97 95 122
+516 4 2 1 1 89 118 87 86
+517 4 2 1 1 15 38 7 81
+518 4 2 1 1 102 44 59 9
+519 4 2 1 1 59 44 102 54
+520 4 2 1 1 110 71 69 75
+521 4 2 1 1 79 77 114 83
+522 4 2 1 1 76 69 109 72
+523 4 2 1 1 113 77 80 84
+524 4 2 1 1 93 99 121 95
+525 4 2 1 1 87 117 91 85
+526 4 2 1 1 90 119 43 11
+527 4 2 1 1 80 78 115 82
+528 4 2 1 1 88 119 85 90
+529 4 2 1 1 90 43 119 85
+530 4 2 1 1 46 108 66 61
+531 4 2 1 1 46 66 108 14
+532 4 2 1 1 26 14 6 73
+533 4 2 1 1 44 101 57 12
+534 4 2 1 1 44 57 101 54
+535 4 2 1 1 46 107 68 15
+536 4 2 1 1 46 68 107 61
+537 4 2 1 1 9 33 1 76
+538 4 2 1 1 100 124 48 94
+539 4 2 1 1 100 48 124 16
+540 4 2 1 1 25 70 13 45
+541 4 2 1 1 56 53 104 58
+542 4 2 1 1 61 64 108 66
+543 4 2 1 1 106 62 67 64
+544 4 2 1 1 102 59 54 56
+545 4 2 1 1 32 13 5 97
+546 4 2 1 1 30 16 8 89
+547 4 2 1 1 13 70 25 74
+548 4 2 1 1 97 122 48 13
+549 4 2 1 1 75 110 41 69
+550 4 2 1 1 54 57 101 55
+551 4 2 1 1 97 48 122 94
+552 4 2 1 1 47 118 89 86
+553 4 2 1 1 47 89 118 16
+554 4 2 1 1 17 9 1 76
+555 4 2 1 1 11 21 3 60
+556 4 2 1 1 10 19 2 58
+557 4 2 1 1 75 41 110 10
+558 4 2 1 1 61 63 68 107
+559 4 2 1 1 70 111 45 112
+560 4 2 1 1 78 115 46 116
+561 4 2 1 1 77 114 42 113
+562 4 2 1 1 69 41 109 110
+563 4 2 1 1 107 61 108 46
+564 4 2 1 1 103 42 104 53
+565 4 2 1 1 118 86 47 120
+566 4 2 1 1 94 124 48 122
+567 4 2 1 1 117 119 43 85
+568 4 2 1 1 123 44 121 93
+569 4 2 1 1 101 54 102 44
+570 4 2 1 1 48 105 62 106
+571 4 2 1 1 91 43 12 22
+572 4 2 1 1 91 43 22 85
+573 4 2 1 1 93 12 23 98
+574 4 2 1 1 93 23 12 44
+575 4 2 1 1 46 81 116 15
+576 4 2 1 1 46 116 81 78
+577 4 2 1 1 92 47 120 15
+578 4 2 1 1 120 47 92 86
+579 4 2 1 1 46 82 115 78
+580 4 2 1 1 46 115 82 14
+581 4 2 1 1 73 45 112 70
+582 4 2 1 1 73 112 45 14
+583 4 2 1 1 99 44 9 24
+584 4 2 1 1 99 44 24 93
+585 4 2 1 1 84 42 113 77
+586 4 2 1 1 84 113 42 10
+587 4 2 1 1 76 41 109 69
+588 4 2 1 1 109 41 76 9
+589 4 2 1 1 42 83 114 77
+590 4 2 1 1 42 114 83 11
+591 4 2 2 2 135 13 122 131
+592 4 2 2 2 132 138 125 110
+593 4 2 2 2 138 108 125 106
+594 4 2 2 2 136 108 125 137
+595 4 2 2 2 110 138 125 131
+596 4 2 2 2 136 107 134 125
+597 4 2 2 2 122 135 131 125
+598 4 2 2 2 13 49 122 131
+599 4 2 2 2 136 108 137 46
+600 4 2 2 2 110 138 131 111
+601 4 2 2 2 132 138 110 112
+602 4 2 2 2 106 107 108 125
+603 4 2 2 2 138 112 137 45
+604 4 2 2 2 104 127 125 128
+605 4 2 2 2 113 115 114 125
+606 4 2 2 2 131 122 125 121
+607 4 2 2 2 112 45 14 137
+608 4 2 2 2 110 131 125 130
+609 4 2 2 2 115 116 114 125
+610 4 2 2 2 132 138 112 137
+611 4 2 2 2 50 137 115 132
+612 4 2 2 2 104 127 128 42
+613 4 2 2 2 105 134 125 133
+614 4 2 2 2 129 120 51 119
+615 4 2 2 2 104 130 128 125
+616 4 2 2 2 105 134 133 47
+617 4 2 2 2 135 106 138 125
+618 4 2 2 2 106 105 107 125
+619 4 2 2 2 138 108 106 45
+620 4 2 2 2 104 130 125 102
+621 4 2 2 2 104 127 42 103
+622 4 2 2 2 109 131 49 111
+623 4 2 2 2 106 135 48 105
+624 4 2 2 2 126 101 43 125
+625 4 2 2 2 138 137 125 108
+626 4 2 2 2 104 127 103 125
+627 4 2 2 2 116 114 129 51
+628 4 2 2 2 136 108 46 107
+629 4 2 2 2 136 108 107 125
+630 4 2 2 2 52 124 123 125
+631 4 2 2 2 132 112 110 50
+632 4 2 2 2 44 121 125 123
+633 4 2 2 2 131 122 121 49
+634 4 2 2 2 105 134 47 107
+635 4 2 2 2 106 13 135 138
+636 4 2 2 2 52 133 125 118
+637 4 2 2 2 103 43 101 125
+638 4 2 2 2 117 43 119 125
+639 4 2 2 2 105 134 107 125
+640 4 2 2 2 137 115 132 125
+641 4 2 2 2 50 137 132 112
+642 4 2 2 2 106 135 105 125
+643 4 2 2 2 102 101 44 125
+644 4 2 2 2 104 130 102 41
+645 4 2 2 2 138 45 13 111
+646 4 2 2 2 131 138 135 13
+647 4 2 2 2 138 112 45 111
+648 4 2 2 2 131 138 13 111
+649 4 2 2 2 107 134 15 136
+650 4 2 2 2 118 52 117 125
+651 4 2 2 2 101 126 12 44
+652 4 2 2 2 104 130 41 128
+653 4 2 2 2 138 137 108 45
+654 4 2 2 2 133 118 16 47
+655 4 2 2 2 44 101 126 125
+656 4 2 2 2 133 134 118 47
+657 4 2 2 2 133 134 125 118
+658 4 2 2 2 124 122 123 125
+659 4 2 2 2 129 118 119 125
+660 4 2 2 2 102 104 103 125
+661 4 2 2 2 135 48 13 106
+662 4 2 2 2 123 122 121 125
+663 4 2 2 2 103 127 43 125
+664 4 2 2 2 106 45 13 138
+665 4 2 2 2 103 11 127 42
+666 4 2 2 2 127 51 129 114
+667 4 2 2 2 13 49 131 111
+668 4 2 2 2 137 108 14 46
+669 4 2 2 2 15 47 134 107
+670 4 2 2 2 10 41 128 104
+671 4 2 2 2 16 118 133 52
+672 4 2 2 2 46 107 15 136
+673 4 2 2 2 41 102 9 130
+674 4 2 2 2 128 50 132 110
+675 4 2 2 2 10 41 110 128
+676 4 2 2 2 50 137 112 14
+677 4 2 2 2 130 102 9 44
+678 4 2 2 2 105 133 16 47
+679 4 2 2 2 127 11 103 43
+680 4 2 2 2 128 130 41 110
+681 4 2 2 2 116 134 15 51
+682 4 2 2 2 137 45 14 108
+683 4 2 2 2 12 126 101 43
+684 4 2 2 2 133 48 135 105
+685 4 2 2 2 128 42 10 104
+686 4 2 2 2 131 9 109 49
+687 4 2 2 2 118 117 119 125
+688 4 2 2 2 102 103 101 125
+689 4 2 2 2 129 118 125 134
+690 4 2 2 2 125 129 116 114
+691 4 2 2 2 117 126 43 125
+692 4 2 2 2 126 117 12 52
+693 4 2 2 2 52 126 117 125
+694 4 2 2 2 12 117 126 43
+695 4 2 2 2 127 51 114 11
+696 4 2 2 2 114 129 127 125
+697 4 2 2 2 127 119 43 125
+698 4 2 2 2 119 11 127 43
+699 4 2 2 2 127 113 114 125
+700 4 2 2 2 127 113 42 114
+701 4 2 2 2 128 130 110 125
+702 4 2 2 2 15 47 120 134
+703 4 2 2 2 127 11 114 42
+704 4 2 2 2 120 47 118 134
+705 4 2 2 2 44 130 102 125
+706 4 2 2 2 44 126 12 123
+707 4 2 2 2 123 44 126 125
+708 4 2 2 2 129 118 134 120
+709 4 2 2 2 110 132 128 125
+710 4 2 2 2 128 127 125 113
+711 4 2 2 2 128 127 113 42
+712 4 2 2 2 128 50 110 10
+713 4 2 2 2 113 42 10 128
+714 4 2 2 2 129 134 125 116
+715 4 2 2 2 127 51 11 119
+716 4 2 2 2 129 134 116 51
+717 4 2 2 2 127 51 119 129
+718 4 2 2 2 129 119 127 125
+719 4 2 2 2 110 131 130 109
+720 4 2 2 2 109 130 9 131
+721 4 2 2 2 121 130 9 44
+722 4 2 2 2 122 48 13 135
+723 4 2 2 2 44 121 130 125
+724 4 2 2 2 110 138 111 112
+725 4 2 2 2 131 138 125 135
+726 4 2 2 2 115 137 14 46
+727 4 2 2 2 126 52 12 123
+728 4 2 2 2 123 126 52 125
+729 4 2 2 2 50 137 14 115
+730 4 2 2 2 136 137 115 46
+731 4 2 2 2 121 9 131 49
+732 4 2 2 2 121 131 130 125
+733 4 2 2 2 131 130 9 121
+734 4 2 2 2 113 132 115 125
+735 4 2 2 2 113 50 115 132
+736 4 2 2 2 128 50 10 113
+737 4 2 2 2 132 113 128 125
+738 4 2 2 2 128 50 113 132
+739 4 2 2 2 133 105 135 125
+740 4 2 2 2 52 124 125 133
+741 4 2 2 2 133 48 105 16
+742 4 2 2 2 16 133 124 52
+743 4 2 2 2 136 137 125 115
+744 4 2 2 2 132 138 137 125
+745 4 2 2 2 41 130 9 109
+746 4 2 2 2 134 120 15 51
+747 4 2 2 2 129 118 120 119
+748 4 2 2 2 129 120 134 51
+749 4 2 2 2 135 122 124 125
+750 4 2 2 2 135 48 124 122
+751 4 2 2 2 133 48 16 124
+752 4 2 2 2 133 135 124 125
+753 4 2 2 2 133 48 124 135
+754 4 2 2 2 116 136 134 125
+755 4 2 2 2 115 136 116 125
+756 4 2 2 2 46 115 136 116
+757 4 2 2 2 136 134 15 116
+758 4 2 2 2 46 136 15 116
+759 4 2 2 2 109 41 130 110
+760 4 2 2 2 110 131 109 111
+$EndElements
@@ -0,0 +1,77 @@
+// Square-in-square 2D geometry for MFEM
+// Creates concentric squares with different material attributes
+
+// Define the square sizes
+L_outer = 2.0;
+L_inner = 0.5;
+
+// Set mesh size and algorithm
+mesh_size = 1.0;
+Mesh.Algorithm = 6; // Frontal-Delaunay for 2D triangular mesh
+Mesh.CharacteristicLengthFactor = 1.0;
+Mesh.MshFileVersion = 2.2;
+
+// Define center point for concentric squares
+cx = 0.0;
+cy = 0.0;
+
+// Define the points (vertices of the outer square)
+Point(1) = {cx-L_outer/2, cy-L_outer/2, 0, mesh_size}; // bottom-left outer
+Point(2) = {cx+L_outer/2, cy-L_outer/2, 0, mesh_size}; // bottom-right outer
+Point(3) = {cx+L_outer/2, cy+L_outer/2, 0, mesh_size}; // top-right outer
+Point(4) = {cx-L_outer/2, cy+L_outer/2, 0, mesh_size}; // top-left outer
+
+// Define the points (vertices of the inner square)
+Point(5) = {cx-L_inner/2, cy-L_inner/2, 0, mesh_size}; // bottom-left inner
+Point(6) = {cx+L_inner/2, cy-L_inner/2, 0, mesh_size}; // bottom-right inner
+Point(7) = {cx+L_inner/2, cy+L_inner/2, 0, mesh_size}; // top-right inner
+Point(8) = {cx-L_inner/2, cy+L_inner/2, 0, mesh_size}; // top-left inner
+
+// Define the lines (edges of the outer square)
+Line(1) = {1, 2}; // bottom edge
+Line(2) = {2, 3}; // right edge
+Line(3) = {3, 4}; // top edge
+Line(4) = {4, 1}; // left edge
+
+// Define the lines (edges of the inner square)
+Line(5) = {5, 6}; // bottom edge
+Line(6) = {6, 7}; // right edge
+Line(7) = {7, 8}; // top edge
+Line(8) = {8, 5}; // left edge
+
+// Define the surfaces
+// Outer square boundary
+Line Loop(1) = {1, 2, 3, 4};
+
+// Inner square boundary (hole in the outer region)
+Line Loop(2) = {5, 6, 7, 8};
+
+// Define the surface areas
+// Outer region (annular region between squares)
+Plane Surface(1) = {1, 2}; // Outer loop minus inner loop (creates hole)
+
+// Inner region (solid inner square)
+Plane Surface(2) = {2}; // Inner loop only
+
+// Assign physical groups for materials
+Physical Surface(1) = {1}; // Outer material (annular region)
+Physical Surface(2) = {2}; // Inner material (solid square)
+
+// Physical lines for boundary conditions
+// Outer square boundary edges
+Physical Line(1) = {1}; // outer bottom
+Physical Line(2) = {2}; // outer right
+Physical Line(3) = {3}; // outer top
+Physical Line(4) = {4}; // outer left
+
+// Inner square boundary edges
+Physical Line(5) = {5}; // inner bottom
+Physical Line(6) = {6}; // inner right
+Physical Line(7) = {7}; // inner top
+Physical Line(8) = {8}; // inner left
+
+// Mesh control for quality
+Mesh.OptimizeNetgen = 1;
+Mesh.Optimize = 1;
+Mesh.ElementOrder = 1;
+Mesh.RecombineAll = 0; // Keep triangular elements (don't recombine to quads)
@@ -0,0 +1,50 @@
+$MeshFormat
+2.2 0 8
+$EndMeshFormat
+$Nodes
+13
+1 -1 -1 0
+2 1 -1 0
+3 1 1 0
+4 -1 1 0
+5 -0.25 -0.25 0
+6 0.25 -0.25 0
+7 0.25 0.25 0
+8 -0.25 0.25 0
+9 -2.752797989558076e-12 -1 0
+10 1 -2.752797989558076e-12 0
+11 2.752797989558076e-12 1 0
+12 -1 2.752797989558076e-12 0
+13 0 0 0
+$EndNodes
+$Elements
+28
+1 1 2 1 1 1 9
+2 1 2 1 1 9 2
+3 1 2 2 2 2 10
+4 1 2 2 2 10 3
+5 1 2 3 3 3 11
+6 1 2 3 3 11 4
+7 1 2 4 4 4 12
+8 1 2 4 4 12 1
+9 1 2 5 5 5 6
+10 1 2 6 6 6 7
+11 1 2 7 7 7 8
+12 1 2 8 8 8 5
+13 2 2 1 1 6 5 9
+14 2 2 1 1 5 8 12
+15 2 2 1 1 7 6 10
+16 2 2 1 1 8 7 11
+17 2 2 1 1 9 5 1
+18 2 2 1 1 5 12 1
+19 2 2 1 1 6 9 2
+20 2 2 1 1 10 6 2
+21 2 2 1 1 7 10 3
+22 2 2 1 1 11 7 3
+23 2 2 1 1 8 11 4
+24 2 2 1 1 8 4 12
+25 2 2 2 2 5 6 13
+26 2 2 2 2 8 5 13
+27 2 2 2 2 6 7 13
+28 2 2 2 2 7 8 13
+$EndElements
@@ -230,11 +230,6 @@ if (MFEM_USE_GINKGO)
  add_subdirectory(ginkgo)
 endif()

-# Include the examples/glvis directory if GLVis is enabled.
-if (MFEM_USE_GLVIS)
-  add_subdirectory(glvis)
-endif()
-
 # Include the examples/hiop directory if HiOp is enabled
 if (MFEM_USE_HIOP)
  add_subdirectory(hiop)
@@ -1,60 +0,0 @@
-# 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.
-
-set(GLVIS_EXAMPLES_SRCS ex0.cpp)
-if (MFEM_USE_MPI)
-  list(APPEND GLVIS_EXAMPLES_SRCS ex0p.cpp)
-endif()
-
-# Include the source directory where mfem.hpp and mfem-performance.hpp are.
-include_directories(BEFORE ${PROJECT_BINARY_DIR})
-
-# Add "test_glvis" target, see below.
-add_custom_target(test_glvis
-  ${CMAKE_CTEST_COMMAND} -R glvis USES_TERMINAL)
-
-# Add one executable per cpp file, adding "glvis_" as prefix.
-# Sets "test_glvis" as a target that depends on the given examples.
-set(PFX glvis_)
-add_mfem_examples(GLVIS_EXAMPLES_SRCS ${PFX} "" test_glvis)
-
-# Testing.
-# The GLVis tests can be run separately using the target "test_glvis"
-# which builds the examples and runs:
-#   ctest -R glvis
-
-if (MFEM_ENABLE_TESTING)
-  # Command line options for the tests.
-  set(EX0_TEST_OPTS -m ../../data/square-disc.mesh)
-  set(EX0P_TEST_OPTS ${EX0_TEST_OPTS})
-
-  # Add the tests: one test per source file.
-  foreach(SRC_FILE ${GLVIS_EXAMPLES_SRCS})
-    get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
-    string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})
-    string(TOUPPER ${TEST_NAME} UP_TEST_NAME)
-    set(TEST_NAME ${PFX}${TEST_NAME})
-
-    set(THIS_TEST_OPTIONS ${${UP_TEST_NAME}_TEST_OPTS})
-    # message(STATUS "Test ${TEST_NAME} options: ${THIS_TEST_OPTIONS}")
-
-    if (NOT (${TEST_NAME} MATCHES ".*p$"))
-      add_test(NAME ${TEST_NAME}_ser
-        COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
-    else()
-      add_test(NAME ${TEST_NAME}_np=${MFEM_MPI_NP}
-        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-        ${MPIEXEC_PREFLAGS}
-        $<TARGET_FILE:${TEST_NAME}> ${THIS_TEST_OPTIONS}
-        ${MPIEXEC_POSTFLAGS})
-    endif()
-  endforeach()
-endif()
@@ -1,27 +0,0 @@
-                    Finite Element Discretization Library
-                                   __
-                       _ __ ___   / _|  ___  _ __ ___
-                      | '_ ` _ \ | |_  / _ \| '_ ` _ \
-                      | | | | | ||  _||  __/| | | | | |
-                      |_| |_| |_||_|   \___||_| |_| |_|
-
-                               https://mfem.org
-
-This directory contains modifications of the example codes that illustrate the
-use of MFEM features based on GLVis in non-server mode.
-
-To build these examples, make sure that MFEM is configured with the option
-"MFEM_USE_GLVIS = YES", see the top-level INSTALL file for details.
-
-Unlike the main examples in examples/, which save files or send data to a
-GLVis server, the codes here use the mfem::glvis_stream class to stream mesh
-and solution data directly into GLVis.
-Currently this directory contains serial and parallel versions of example 0
-(ex0 and ex0p).
-
-We recommend comparing the original example codes with the corresponding files
-in the current directory.
-
-From this directory, the codes can be built with "make" and tested with
-"make test". With CMake, the executables are named glvis_ex0 and glvis_ex0p,
-and the tests can be run with "make test_glvis" or "ctest -R glvis".
@@ -1,82 +0,0 @@
-//                                MFEM Example 0
-//
-// Compile with: make ex0
-//
-// Sample runs:  ex0
-//               ex0 -m ../data/fichera.mesh
-//               ex0 -m ../data/square-disc.mesh -o 2
-//
-// Description: This example code demonstrates the most basic usage of MFEM to
-//              define a simple finite element discretization of the Poisson
-//              problem -Delta u = 1 with zero Dirichlet boundary conditions.
-//              General 2D/3D mesh files and finite element polynomial degrees
-//              can be specified by command line options.
-
-#include "mfem.hpp"
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command line options.
-   string mesh_file = "../data/star.mesh";
-   int order = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
-   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
-   args.ParseCheck();
-
-   // 2. Read the mesh from the given mesh file, and refine once uniformly.
-   Mesh mesh(mesh_file);
-   mesh.UniformRefinement();
-
-   // 3. Define a finite element space on the mesh. Here we use H1 continuous
-   //    high-order Lagrange finite elements of the given order.
-   H1_FECollection fec(order, mesh.Dimension());
-   FiniteElementSpace fespace(&mesh, &fec);
-   cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;
-
-   // 4. Extract the list of all the boundary DOFs. These will be marked as
-   //    Dirichlet in order to enforce zero boundary conditions.
-   Array<int> boundary_dofs;
-   fespace.GetBoundaryTrueDofs(boundary_dofs);
-
-   // 5. Define the solution x as a finite element grid function in fespace. Set
-   //    the initial guess to zero, which also sets the boundary conditions.
-   GridFunction x(&fespace);
-   x = 0.0;
-
-   // 6. Set up the linear form b(.) corresponding to the right-hand side.
-   ConstantCoefficient one(1.0);
-   LinearForm b(&fespace);
-   b.AddDomainIntegrator(new DomainLFIntegrator(one));
-   b.Assemble();
-
-   // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
-   BilinearForm a(&fespace);
-   a.AddDomainIntegrator(new DiffusionIntegrator);
-   a.Assemble();
-
-   // 8. Form the linear system A X = B. This includes eliminating boundary
-   //    conditions, applying AMR constraints, and other transformations.
-   SparseMatrix A;
-   Vector B, X;
-   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
-
-   // 9. Solve the system using PCG with symmetric Gauss-Seidel preconditioner.
-   GSSmoother M(A);
-   PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
-
-   // 10. Recover the solution x as a grid function and save to file.
-   a.RecoverFEMSolution(X, b, x);
-
-   // 11. Send the solution to a non-server mode GLVis.
-   glvis_stream glvis;
-   glvis.precision(8);
-   glvis << "solution\n" << mesh << x << flush;
-
-   return 0;
-}
@@ -1,103 +0,0 @@
-//                       MFEM Example 0 - Parallel Version
-//
-// Compile with: make ex0p
-//
-// Sample runs:  mpirun -np 4 ex0p
-//               mpirun -np 4 ex0p -m ../data/fichera.mesh
-//               mpirun -np 4 ex0p -m ../data/square-disc.mesh -o 2
-//
-// Description: This example code demonstrates the most basic parallel usage of
-//              MFEM to define a simple finite element discretization of the
-//              Poisson problem -Delta u = 1 with zero Dirichlet boundary
-//              conditions. General 2D/3D serial mesh files and finite element
-//              polynomial degrees can be specified by command line options.
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI and HYPRE.
-   Mpi::Init(argc, argv);
-   Hypre::Init();
-
-   // 2. Parse command line options.
-   string mesh_file = "../data/star.mesh";
-   int order = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
-   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
-   args.ParseCheck();
-
-   // 3. Read the serial mesh from the given mesh file.
-   Mesh serial_mesh(mesh_file);
-
-   // 4. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh once in parallel to increase the resolution.
-   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
-   serial_mesh.Clear(); // the serial mesh is no longer needed
-   mesh.UniformRefinement();
-
-   // 5. Define a finite element space on the mesh. Here we use H1 continuous
-   //    high-order Lagrange finite elements of the given order.
-   H1_FECollection fec(order, mesh.Dimension());
-   ParFiniteElementSpace fespace(&mesh, &fec);
-   HYPRE_BigInt total_num_dofs = fespace.GlobalTrueVSize();
-   if (Mpi::Root())
-   {
-      cout << "Number of unknowns: " << total_num_dofs << endl;
-   }
-
-   // 6. Extract the list of all the boundary DOFs. These will be marked as
-   //    Dirichlet in order to enforce zero boundary conditions.
-   Array<int> boundary_dofs;
-   fespace.GetBoundaryTrueDofs(boundary_dofs);
-
-   // 7. Define the solution x as a finite element grid function in fespace. Set
-   //    the initial guess to zero, which also sets the boundary conditions.
-   ParGridFunction x(&fespace);
-   x = 0.0;
-
-   // 8. Set up the linear form b(.) corresponding to the right-hand side.
-   ConstantCoefficient one(1.0);
-   ParLinearForm b(&fespace);
-   b.AddDomainIntegrator(new DomainLFIntegrator(one));
-   b.Assemble();
-
-   // 9. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
-   ParBilinearForm a(&fespace);
-   a.AddDomainIntegrator(new DiffusionIntegrator);
-   a.Assemble();
-
-   // 10. Form the linear system A X = B. This includes eliminating boundary
-   //     conditions, applying AMR constraints, parallel assembly, etc.
-   HypreParMatrix A;
-   Vector B, X;
-   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
-
-   // 11. Solve the system using PCG with hypre's BoomerAMG preconditioner.
-   HypreBoomerAMG M(A);
-   CGSolver cg(MPI_COMM_WORLD);
-   cg.SetRelTol(1e-12);
-   cg.SetMaxIter(2000);
-   cg.SetPrintLevel(1);
-   cg.SetPreconditioner(M);
-   cg.SetOperator(A);
-   cg.Mult(B, X);
-
-   // 12. Recover the solution x as a grid function.
-   a.RecoverFEMSolution(X, b, x);
-
-   // 13. Send the solution to a non-server mode GLVis.
-   glvis_stream glvis;
-   glvis << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n";
-   glvis.precision(8);
-   glvis << "solution\n" << mesh << x << flush;
-
-   return 0;
-}
@@ -1,74 +0,0 @@
-# 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.
-
-# Use the MFEM build directory
-MFEM_DIR ?= ../..
-MFEM_BUILD_DIR ?= ../..
-MFEM_INSTALL_DIR ?= ../../mfem
-SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/glvis/,)
-CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
-   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
-
-MFEM_LIB_FILE = mfem_is_not_built
-include $(CONFIG_MK)
-
-SEQ_EXAMPLES = ex0
-PAR_EXAMPLES = ex0p
-ifeq ($(MFEM_USE_MPI),NO)
-   EXAMPLES = $(SEQ_EXAMPLES)
-else
-   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
-endif
-
-.SUFFIXES:
-.SUFFIXES: .o .cpp .mk
-.PHONY: all clean clean-build clean-exec
-
-# Remove built-in rule
-%: %.cpp
-
-# Replace the default implicit rule for *.cpp files
-%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
-
-all: $(EXAMPLES)
-
-ifeq ($(MFEM_USE_GLVIS),NO)
-$(EXAMPLES):
-	$(error MFEM is not configured with GLVIS)
-endif
-
-MFEM_TESTS = EXAMPLES
-include $(MFEM_TEST_MK)
-
-# Testing
-RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
-EX0_ARGS := -m ../../data/square-disc.mesh
-
-ex0-test-seq: ex0
-	@$(call mfem-test,$<,, Serial GLVis example,$(EX0_ARGS),SKIP-NO-VIS)
-ex0p-test-par: ex0p
-	@$(call mfem-test,$<, $(RUN_MPI), Parallel GLVis example,$(EX0_ARGS),SKIP-NO-VIS)
-
-# Testing: "test" target and mfem-test* variables are defined in config/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 *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
-	rm -rf *.dSYM *.TVD.*breakpoints
-
-clean-exec:
-	@:
@@ -22,6 +22,8 @@
 #include <algorithm>
 #include <cmath>
 #include <cstdarg>
+#include <unordered_map>
+#include <unordered_set>

 using namespace std;

@@ -4527,6 +4529,210 @@ void FiniteElementSpace
   }
 }

+void FiniteElementSpace::GetBoundaryLoopEdgeDofs(
+   const Array<int> &boundary_element_indices,
+   Array<int> &boundary_edge_dofs,
+   Array<int> *dof_edges,
+   Array<int> *dof_boundary_elements) const
+{
+   MFEM_VERIFY(mesh->Dimension() >= 2,
+               "GetBoundaryLoopEdgeDofs requires 2D or 3D meshes to find edge objects");
+
+   boundary_edge_dofs.SetSize(0);
+   if (dof_edges) { dof_edges->SetSize(0); }
+   if (dof_boundary_elements) { dof_boundary_elements->SetSize(0); }
+
+   // A DOF that appears in exactly one selected boundary element lies on the
+   // bounding loop; one appearing in two or more is interior to the boundary
+   // region and is dropped. Count occurrences of each DOF (using scratch maps,
+   // exposed only as parallel-indexed Array<int> below) and record, on first
+   // sight, the local edge and boundary element carrying it.
+   //
+   // The count is over GetEdgeDofs, which returns endpoint vertex DOFs as well
+   // as edge-interior DOFs (relevant for collections such as ND_R2D that carry
+   // vertex DOFs). Edge-interior DOFs occur once per edge, so the count mainly
+   // resolves vertex DOFs: a vertex shared by several elements is interior and
+   // dropped, while a genuine loop-corner (open-curve endpoint) vertex is kept.
+   // This is why we count GetEdgeDofs rather than collecting GetEdgeInteriorDofs,
+   // which would omit the endpoint vertex DOFs the method is documented to keep.
+   // The 3D removal criterion (any edge in two or more faces) matches the
+   // parallel version rather than a parity toggle.
+   std::unordered_map<int, int> dof_count, dof_edge, dof_belem;
+   Array<int> edge_dofs, edges, edge_orientations;
+
+   const int dim = mesh->Dimension();
+   for (int i = 0; i < boundary_element_indices.Size(); ++i)
+   {
+      const int boundary_element_idx = boundary_element_indices[i];
+      std::unordered_set<int> boundary_element_dofs;
+
+      if (dim == 3)
+      {
+         // Boundary elements are 2D faces; extract their 1D edges.
+         int face_index, face_orientation;
+         mesh->GetBdrElementFace(boundary_element_idx, &face_index,
+                                 &face_orientation);
+         mesh->GetFaceEdges(face_index, edges, edge_orientations);
+      }
+      else
+      {
+         // Boundary elements are 1D segments, each being a single edge.
+         mesh->GetBdrElementEdges(boundary_element_idx, edges, edge_orientations);
+         MFEM_VERIFY(edges.Size() == 1,
+                     "2D boundary element should have exactly one edge");
+      }
+
+      for (int j = 0; j < edges.Size(); ++j)
+      {
+         GetEdgeDofs(edges[j], edge_dofs);
+         for (int k = 0; k < edge_dofs.Size(); ++k)
+         {
+            const int dof = edge_dofs[k];
+            // Count each DOF once per boundary element and record metadata the
+            // first time it is seen, so H1 DOFs shared by multiple edges of the
+            // same element are not double counted.
+            if (boundary_element_dofs.insert(dof).second &&
+                dof_count[dof]++ == 0)
+            {
+               dof_edge[dof] = edges[j];
+               dof_belem[dof] = boundary_element_idx;
+            }
+         }
+      }
+   }
+
+   // Emit the DOFs seen in exactly one selected boundary element, in a
+   // deterministic (increasing DOF index) order shared by all output arrays.
+   std::vector<int> kept;
+   kept.reserve(dof_count.size());
+   for (const auto &[dof, count] : dof_count)
+   {
+      if (count == 1) { kept.push_back(dof); }
+   }
+   std::sort(kept.begin(), kept.end());
+
+   boundary_edge_dofs.Reserve(static_cast<int>(kept.size()));
+   if (dof_edges) { dof_edges->Reserve(static_cast<int>(kept.size())); }
+   if (dof_boundary_elements)
+   {
+      dof_boundary_elements->Reserve(static_cast<int>(kept.size()));
+   }
+   for (int dof : kept)
+   {
+      boundary_edge_dofs.Append(dof);
+      if (dof_edges) { dof_edges->Append(dof_edge[dof]); }
+      if (dof_boundary_elements) { dof_boundary_elements->Append(dof_belem[dof]); }
+   }
+}
+
+void FiniteElementSpace::GetBoundaryElementsByAttribute(
+   const Array<int> &bdr_attrs,
+   std::vector<Array<int>> &attr_to_elements)
+{
+   // One (initially empty) list of boundary elements per requested attribute,
+   // indexed to match bdr_attrs.
+   attr_to_elements.assign(bdr_attrs.Size(), Array<int>());
+
+   // Map attribute value -> position in bdr_attrs for quick lookup.
+   std::unordered_map<int, int> attr_to_index;
+   for (int i = 0; i < bdr_attrs.Size(); ++i)
+   {
+      attr_to_index[bdr_attrs[i]] = i;
+   }
+
+   // Bucket boundary elements by their attribute.
+   for (int i = 0; i < mesh->GetNBE(); ++i)
+   {
+      int attr = mesh->GetBdrElement(i)->GetAttribute();
+      auto it = attr_to_index.find(attr);
+      if (it != attr_to_index.end())
+      {
+         attr_to_elements[it->second].Append(i);
+      }
+   }
+}
+
+void FiniteElementSpace::GetBoundaryElementsByAttribute(int bdr_attr,
+                                                        Array<int> &boundary_elements)
+{
+   boundary_elements.SetSize(0);
+
+   for (int i = 0; i < mesh->GetNBE(); ++i)
+   {
+      if (mesh->GetBdrElement(i)->GetAttribute() == bdr_attr)
+      {
+         boundary_elements.Append(i);
+      }
+   }
+}
+
+void FiniteElementSpace::ComputeLoopEdgeOrientations(
+   const Array<int> &dof_edges,
+   const Array<int> &dof_boundary_elements,
+   const Vector &loop_normal,
+   Array<int> &dof_orientations) const
+{
+   MFEM_VERIFY(dof_edges.Size() == dof_boundary_elements.Size(),
+               "dof_edges and dof_boundary_elements must be parallel-indexed");
+
+   const int ndof = dof_edges.Size();
+   dof_orientations.SetSize(ndof);
+
+   Array<int> edge_verts, bdr_elem_verts;
+   Vector edge_vec(3), to_edge_vec(3), cross_product(3);
+   for (int i = 0; i < ndof; i++)
+   {
+      const int edge_id = dof_edges[i];
+      const int bdr_elem_idx = dof_boundary_elements[i];
+
+      // Get edge vertices
+      mesh->GetEdgeVertices(edge_id, edge_verts);
+
+      const real_t *v0 = mesh->GetVertex(edge_verts[0]);
+      const real_t *v1 = mesh->GetVertex(edge_verts[1]);
+
+      // Get boundary element vertices
+      mesh->GetBdrElement(bdr_elem_idx)->GetVertices(bdr_elem_verts);
+
+      // Find the third vertex (not part of the edge)
+      int third_vertex = -1;
+      for (int j = 0; j < bdr_elem_verts.Size(); j++)
+      {
+         int v = bdr_elem_verts[j];
+         if (v != edge_verts[0] && v != edge_verts[1])
+         {
+            third_vertex = v;
+            break;
+         }
+      }
+
+      if (third_vertex == -1)
+      {
+         MFEM_ABORT("Boundary element " << bdr_elem_idx << " has only 2 vertices, "
+                    "but 3D boundary elements must have at least 3 vertices");
+      }
+
+      const real_t *v2 = mesh->GetVertex(third_vertex);
+
+      // Edge vector
+      for (int j = 0; j < 3; j++) { edge_vec[j] = v1[j] - v0[j]; }
+
+      // Vector from third vertex to edge (use edge midpoint)
+      for (int j = 0; j < 3; j++)
+      {
+         real_t edge_midpoint = (v0[j] + v1[j]) * 0.5;
+         to_edge_vec[j] = edge_midpoint - v2[j];
+      }
+
+      // Cross product: to_edge × edge
+      to_edge_vec.cross3D(edge_vec, cross_product);
+
+      // Check alignment with loop normal
+      real_t dot_product = cross_product * loop_normal;
+      dof_orientations[i] = (dot_product > 0) ? 1 : -1;
+   }
+}
+
 FiniteElementCollection *FiniteElementSpace::Load(Mesh *m, std::istream &input)
 {
   string buff;
@@ -22,6 +22,7 @@
 #include "restriction.hpp"
 #include <iostream>
 #include <unordered_map>
+#include <vector>

 namespace mfem
 {
@@ -1389,6 +1390,80 @@ public:
   virtual void GetExteriorTrueDofs(Array<int> &exterior_dofs,
                                    int component = -1) const;

+   /** @brief Extract the edge degrees of freedom of a boundary "loop".
+
+       Here a "loop" is the set of boundary edges bounding the region covered by
+       @a boundary_element_indices: in 3D the outer edges of a patch of boundary
+       faces, in 2D the boundary segments themselves. An edge that is shared by
+       two (or more) of the selected boundary elements is interior to that region
+       rather than on its bounding loop, so its DOFs are excluded from the result.
+       This exclusion of interior DOFs is the defining feature of the method.
+
+       The three output arrays share a single indexing: for each valid index @a i,
+       @a dof_edges[i] and @a dof_boundary_elements[i] describe the DOF
+       @a boundary_edge_dofs[i].
+
+       @param[in]  boundary_element_indices Boundary element indices spanning a
+                   boundary surface (3D) or curve (2D).
+       @param[out] boundary_edge_dofs Local DOF indices on the boundary loop.
+       @param[out] dof_edges Optional; local edge index carrying each DOF.
+       @param[out] dof_boundary_elements Optional; a boundary element containing
+                   each DOF.
+
+       @note In 3D the edge DOFs are extracted from the 1D edges of the 2D
+       boundary faces; in 2D they come directly from the 1D boundary segments, so
+       @a dof_edges then holds the boundary element (segment) edge indices.
+       @note This method uses GetEdgeDofs internally, which returns both vertex and
+       edge DOFs. Standard Nédélec elements (ND_FECollection) have no vertex DOFs,
+       so only genuine edge DOFs appear. Collections that carry vertex DOFs (e.g.
+       ND_R2D_FECollection) additionally contribute the vertex DOFs at loop
+       endpoints.
+       @note This is the serial version. For parallel meshes, use the parallel
+       version in ParFiniteElementSpace which handles processor boundaries
+       correctly.
+       @note Requires a 2D or 3D mesh to identify edge objects. The method will
+       assert if called on 1D meshes.
+       @note Only supports conforming meshes; non-conforming meshes are not
+       supported. */
+   void GetBoundaryLoopEdgeDofs(const Array<int> &boundary_element_indices,
+                                Array<int> &boundary_edge_dofs,
+                                Array<int> *dof_edges = nullptr,
+                                Array<int> *dof_boundary_elements = nullptr) const;
+
+   /** @brief Get boundary elements grouped by attribute.
+
+       For each attribute in @a bdr_attrs, collect the indices of all boundary
+       elements carrying that attribute. The result is indexed to match
+       @a bdr_attrs: @a attr_to_elements[i] holds the boundary elements with
+       attribute @a bdr_attrs[i]. */
+   void GetBoundaryElementsByAttribute(
+      const Array<int> &bdr_attrs,
+      std::vector<Array<int>> &attr_to_elements);
+
+   /** @brief Get all boundary elements with a specific attribute. */
+   void GetBoundaryElementsByAttribute(int bdr_attr,
+                                       Array<int> &boundary_elements);
+
+   /** @brief Compute edge orientations relative to a boundary loop direction.
+
+       For each boundary-loop DOF described by @a dof_edges and
+       @a dof_boundary_elements (see GetBoundaryLoopEdgeDofs), determine whether
+       the carrying edge is
+       traversed in the direction consistent with @a loop_normal, following the
+       right-hand rule. Intended for 3D meshes.
+
+       @param[in]  dof_edges Local edge index of each DOF (parallel-indexed with
+                   the boundary_edge_dofs output of GetBoundaryLoopEdgeDofs).
+       @param[in]  dof_boundary_elements A boundary element containing each DOF,
+                   using the same indexing as @a dof_edges.
+       @param[in]  loop_normal Normal vector defining the loop orientation.
+       @param[out] dof_orientations Orientation (+1 or -1) for each DOF, using the
+                   same indexing as @a dof_edges. */
+   void ComputeLoopEdgeOrientations(const Array<int> &dof_edges,
+                                    const Array<int> &dof_boundary_elements,
+                                    const Vector &loop_normal,
+                                    Array<int> &dof_orientations) const;
+
   /// Convert a Boolean marker array to a list containing all marked indices.
   static void MarkerToList(const Array<int> &marker, Array<int> &list);

@@ -26,6 +26,8 @@

 #include <limits>
 #include <list>
+#include <unordered_map>
+#include <unordered_set>

 namespace mfem
 {
@@ -1285,6 +1287,342 @@ void ParFiniteElementSpace::GetExteriorVDofs(Array<int> &ext_dofs,
   Synchronize(ext_dofs);
 }

+void ParFiniteElementSpace::GetBoundaryLoopEdgeDofs(
+   const Array<int> &boundary_element_indices,
+   Array<int> &ess_tdof_list,
+   Array<int> &boundary_edge_dofs_out,
+   Array<int> *ldof_marker,
+   Array<int> *dof_edges,
+   Array<int> *dof_boundary_elements,
+   Array<int> *ess_edge_list)
+{
+   MFEM_VERIFY(!pmesh->Nonconforming(),
+               "GetBoundaryLoopEdgeDofs does not support nonconforming meshes");
+   MFEM_VERIFY(pmesh->Dimension() >= 2,
+               "GetBoundaryLoopEdgeDofs requires 2D or 3D meshes to find 1D edge objects");
+
+   // Call the serial version, then rebuild scratch maps/set from the returned
+   // arrays for the O(1) lookups the parallel reconciliation below needs.
+   Array<int> loc_dofs, loc_edges, loc_belems;
+   FiniteElementSpace::GetBoundaryLoopEdgeDofs(boundary_element_indices, loc_dofs,
+                                               &loc_edges, &loc_belems);
+
+   std::unordered_set<int> boundary_edge_dofs;
+   std::unordered_map<int, int> dof_to_edge_map;
+   std::unordered_map<int, int> dof_to_boundary_element;
+   boundary_edge_dofs.reserve(loc_dofs.Size());
+   dof_to_edge_map.reserve(loc_dofs.Size());
+   dof_to_boundary_element.reserve(loc_dofs.Size());
+   for (int i = 0; i < loc_dofs.Size(); i++)
+   {
+      const int dof = loc_dofs[i];
+      boundary_edge_dofs.insert(dof);
+      dof_to_edge_map[dof] = loc_edges[i];
+      dof_to_boundary_element[dof] = loc_belems[i];
+   }
+
+   // Parallel processing: Build edge sharing lookup table
+   std::unordered_map<int, int> edge_to_group_size;
+   int num_groups = pmesh->GetNGroups();
+
+   int total_shared_edges = 0;
+   for (int group = 1; group < num_groups; group++)
+   {
+      total_shared_edges += pmesh->GroupNEdges(group);
+   }
+   edge_to_group_size.reserve(total_shared_edges);
+
+   for (int group = 1; group < num_groups; group++)
+   {
+      int group_size = pmesh->gtopo.GetGroupSize(group);
+      int num_edges_in_group = pmesh->GroupNEdges(group);
+
+      for (int i = 0; i < num_edges_in_group; i++)
+      {
+         edge_to_group_size.emplace(pmesh->GroupEdge(group, i), group_size);
+      }
+   }
+
+   // Get global indices
+   Array<HYPRE_BigInt> global_edge_indices;
+   pmesh->GetGlobalEdgeIndices(global_edge_indices);
+
+   // Handle dimension-specific boundary element relationships
+   Array<HYPRE_BigInt> global_face_indices;
+   std::unordered_map<int, int> boundary_element_to_companion;
+   std::unordered_set<int> dofs_to_remove;
+
+   const int dim = pmesh->Dimension();
+   if (dim == 3)
+   {
+      // In 3D: boundary elements are faces, we track which face each boundary element is
+      pmesh->GetGlobalFaceIndices(global_face_indices);
+      for (int boundary_element_idx : boundary_element_indices)
+      {
+         int face_index, face_orientation;
+         pmesh->GetBdrElementFace(boundary_element_idx, &face_index, &face_orientation);
+         boundary_element_to_companion[boundary_element_idx] = face_index;
+      }
+
+      std::vector<HYPRE_BigInt> local_data;
+      local_data.reserve(boundary_edge_dofs.size() * 2);
+
+      std::unordered_set<int> processed_edges;
+      processed_edges.reserve(boundary_edge_dofs.size());
+
+      for (const auto& [dof, local_edge] : dof_to_edge_map)
+      {
+         // Skip if already processed this edge
+         if (!processed_edges.insert(local_edge).second) { continue; }
+
+         // Check if edge is shared (fast lookup)
+         auto it = edge_to_group_size.find(local_edge);
+         if (it != edge_to_group_size.end() && it->second > 1)
+         {
+            // Get boundary element and companion index directly from pre-computed map
+            int boundary_element_idx = dof_to_boundary_element[dof];
+            int companion_index = boundary_element_to_companion[boundary_element_idx];
+
+            // Store edge-face pair for 3D artificial boundary detection
+            local_data.push_back(global_edge_indices[local_edge]);
+            local_data.push_back(global_face_indices[companion_index]);
+         }
+      }
+
+      // MPI communication for 3D artificial boundary detection
+      int num_procs = pmesh->GetNRanks();
+      int local_size = local_data.size();
+
+      std::vector<int> mpi_arrays(num_procs * 4);
+      int* all_sizes = mpi_arrays.data();
+      int* displs = all_sizes + num_procs;
+      int* byte_sizes = displs + num_procs;
+      int* byte_displs = byte_sizes + num_procs;
+
+      MPI_Allgather(&local_size, 1, MPI_INT, all_sizes, 1, MPI_INT, pmesh->GetComm());
+
+      int total_size = 0;
+      constexpr int hypre_size = sizeof(HYPRE_BigInt);
+      for (int i = 0; i < num_procs; i++)
+      {
+         displs[i] = total_size;
+         byte_displs[i] = total_size * hypre_size;
+         total_size += all_sizes[i];
+         byte_sizes[i] = all_sizes[i] * hypre_size;
+      }
+
+      if (total_size > 0)
+      {
+         std::vector<HYPRE_BigInt> all_data(total_size);
+         MPI_Allgatherv(local_data.data(), local_size * hypre_size, MPI_BYTE,
+                        all_data.data(), byte_sizes, byte_displs, MPI_BYTE, pmesh->GetComm());
+
+         // Build global-to-local edge mapping
+         std::unordered_map<HYPRE_BigInt, int> global_to_local_edge;
+         global_to_local_edge.reserve(global_edge_indices.Size());
+         for (int i = 0; i < global_edge_indices.Size(); ++i)
+         {
+            global_to_local_edge[global_edge_indices[i]] = i;
+         }
+
+         // Process collected data to find edges in multiple faces (artificial boundaries)
+         std::unordered_map<HYPRE_BigInt, std::unordered_set<HYPRE_BigInt>>edge_to_faces;
+         edge_to_faces.reserve(total_size / 2);
+
+         for (size_t i = 0; i < all_data.size(); i += 2)
+         {
+            edge_to_faces[all_data[i]].insert(all_data[i + 1]);
+         }
+
+         // Mark DOFs from artificial edges for removal
+         dofs_to_remove.reserve(local_data.size() / 4);
+
+         for (size_t i = 0; i < local_data.size(); i += 2)
+         {
+            HYPRE_BigInt global_edge_id = local_data[i];
+
+            // If this edge appears in 2+ distinct faces, it's artificial
+            if (edge_to_faces[global_edge_id].size() >= 2)
+            {
+               int local_edge = global_to_local_edge[global_edge_id];
+               Array<int> local_edge_dofs;
+               GetEdgeDofs(local_edge, local_edge_dofs);
+
+               // Mark boundary DOFs of this edge for removal
+               for (int k = 0; k < local_edge_dofs.Size(); ++k)
+               {
+                  int dof = local_edge_dofs[k];
+                  if (boundary_edge_dofs.count(dof))
+                  {
+                     dofs_to_remove.insert(dof);
+                  }
+               }
+            }
+         }
+      }
+   }
+   else if (dim == 2)
+   {
+      // In 2D the boundary elements are themselves the edges, so there are no
+      // artificial boundary edges to detect. However, for collections with
+      // vertex DOFs (e.g. ND_R2D), a vertex shared by two boundary segments is
+      // interior to the boundary curve and must be dropped. The serial code
+      // does this by erasing a DOF on its second occurrence, which only sees
+      // the occurrences local to this rank. When the two segments meeting at a
+      // vertex live on different ranks, each rank sees a single occurrence and
+      // wrongly keeps the DOF. Reconcile the occurrence parity across each
+      // sharing group: membership in boundary_edge_dofs is the local parity,
+      // and the parities sum (mod 2) to the global occurrence parity.
+      Array<int> boundary_dof_count(GetVSize());
+      boundary_dof_count = 0;
+      for (const int dof : boundary_edge_dofs)
+      {
+         boundary_dof_count[dof] = 1;
+      }
+
+      // implement allreduce(+) as reduce(+) + broadcast
+      gcomm->Reduce<int>(boundary_dof_count, GroupCommunicator::Sum);
+      gcomm->Bcast(boundary_dof_count);
+
+      for (const int dof : boundary_edge_dofs)
+      {
+         if (boundary_dof_count[dof] % 2 == 0)
+         {
+            dofs_to_remove.insert(dof);
+         }
+      }
+   }
+
+   // Remove artificial DOFs
+   for (int dof : dofs_to_remove)
+   {
+      boundary_edge_dofs.erase(dof);
+      dof_to_edge_map.erase(dof);
+      dof_to_boundary_element.erase(dof);
+   }
+
+   // Convert to true DOFs and output
+   ess_tdof_list.SetSize(0);
+   ess_tdof_list.Reserve(boundary_edge_dofs.size());
+   if (ess_edge_list)
+   {
+      // Reset as well, so that it stays in correspondence with ess_tdof_list
+      // when the same output array is reused across calls.
+      ess_edge_list->SetSize(0);
+      ess_edge_list->Reserve(boundary_edge_dofs.size());
+   }
+   // Marker of the boundary edge DOFs. Always computed locally because the
+   // parallel reconciliation below needs it; only copied to the caller's output
+   // if requested (see the ldof_marker parameter).
+   Array<int> local_ldof_marker(GetVSize());
+   local_ldof_marker = 0;
+
+   for (int dof : boundary_edge_dofs)
+   {
+      local_ldof_marker[dof] = 1; // Mark all boundary edge dofs
+   }
+
+   // Make sure that a selected shared DOF is marked on every rank of its
+   // sharing group, including ranks holding none of the selected boundary
+   // elements. Only the group master owns the corresponding true DOF, so
+   // without this the true DOF would be emitted by no rank at all: the
+   // non-master ranks get -1 from GetLocalTDofNumber(), while the master may
+   // not have selected the DOF locally.
+   Synchronize(local_ldof_marker);
+
+   // A DOF marked only through the synchronization above has no local
+   // dof_to_edge_map entry, but the shared edge carrying it is still present in
+   // the local mesh. Build the missing DOF -> edge entries from the shared
+   // edges of the groups, so that ess_edge_list stays in correspondence with
+   // ess_tdof_list. Note that a vertex DOF is not associated with a unique
+   // edge, so it is only resolved when it is an interior DOF of an edge.
+   std::unordered_map<int, int> shared_dof_to_edge;
+   Array<int> shared_edge_dofs;
+   for (int group = 1; group < num_groups; group++)
+   {
+      const int num_edges_in_group = pmesh->GroupNEdges(group);
+      for (int i = 0; i < num_edges_in_group; i++)
+      {
+         const int edge = pmesh->GroupEdge(group, i);
+         GetEdgeInteriorDofs(edge, shared_edge_dofs);
+         for (int k = 0; k < shared_edge_dofs.Size(); k++)
+         {
+            shared_dof_to_edge.emplace(shared_edge_dofs[k], edge);
+         }
+      }
+   }
+
+   // Build parallel arrays for DOFs and corresponding edges
+   std::vector<std::pair<int, int>> tdof_edge_pairs;
+   tdof_edge_pairs.reserve(boundary_edge_dofs.size());
+
+   for (int dof = 0; dof < local_ldof_marker.Size(); dof++)
+   {
+      if (!local_ldof_marker[dof]) { continue; }
+
+      const int tdof = GetLocalTDofNumber(dof);
+      if (tdof < 0) { continue; } // tdof == -1 means not owned by this rank
+
+      int edge = -1;
+      auto it = dof_to_edge_map.find(dof);
+      if (it != dof_to_edge_map.end())
+      {
+         edge = it->second;
+      }
+      else
+      {
+         auto shared_it = shared_dof_to_edge.find(dof);
+         if (shared_it != shared_dof_to_edge.end())
+         {
+            edge = shared_it->second;
+         }
+      }
+      tdof_edge_pairs.push_back({tdof, edge});
+   }
+
+   // Sort by true DOF index to maintain consistent ordering
+   std::sort(tdof_edge_pairs.begin(), tdof_edge_pairs.end());
+
+   // Extract sorted true DOFs and edges
+   for (const auto& pair : tdof_edge_pairs)
+   {
+      ess_tdof_list.Append(pair.first);
+      if (ess_edge_list)
+      {
+         ess_edge_list->Append(pair.second);
+      }
+   }
+
+   // Emit the local boundary-loop DOFs in a deterministic (increasing DOF
+   // index) order shared by all output arrays.
+   std::vector<int> kept(boundary_edge_dofs.begin(), boundary_edge_dofs.end());
+   std::sort(kept.begin(), kept.end());
+
+   boundary_edge_dofs_out.SetSize(0);
+   boundary_edge_dofs_out.Reserve(static_cast<int>(kept.size()));
+   if (dof_edges)
+   {
+      dof_edges->SetSize(0);
+      dof_edges->Reserve(static_cast<int>(kept.size()));
+   }
+   if (dof_boundary_elements)
+   {
+      dof_boundary_elements->SetSize(0);
+      dof_boundary_elements->Reserve(static_cast<int>(kept.size()));
+   }
+   for (int dof : kept)
+   {
+      boundary_edge_dofs_out.Append(dof);
+      if (dof_edges) { dof_edges->Append(dof_to_edge_map[dof]); }
+      if (dof_boundary_elements)
+      {
+         dof_boundary_elements->Append(dof_to_boundary_element[dof]);
+      }
+   }
+
+   if (ldof_marker) { ldof_marker->Swap(local_ldof_marker); }
+}
+
 void ParFiniteElementSpace::GetExteriorTrueDofs(Array<int> &ext_tdof_list,
                                                int component) const
 {
@@ -460,6 +460,41 @@ public:
   void GetExteriorTrueDofs(Array<int> &ext_tdof_list,
                            int component = -1) const override;

+   /** @brief Extract the edge degrees of freedom of a boundary "loop" on a
+       parallel mesh (see the serial FiniteElementSpace::GetBoundaryLoopEdgeDofs
+       for the definition of a loop). This version removes the artificial
+       boundary edges that appear at processor boundaries, so the selected DOFs
+       are independent of the mesh partitioning.
+
+       As in the serial version, the @a boundary_edge_dofs_out, @a dof_edges and
+       @a dof_boundary_elements outputs share a single indexing describing the
+       same local DOF at each position.
+
+       Requirements:
+       - Mesh must be conforming (no hanging nodes)
+       - Mesh dimension must be >= 2
+       @param[in]  boundary_element_indices Array of boundary element indices.
+       @param[out] ess_tdof_list Essential true DOF indices, sorted ascending.
+       @param[out] boundary_edge_dofs_out Local boundary-loop DOF indices.
+       @param[out] ldof_marker Optional; marker of the boundary edge DOFs,
+                   derivable from @a boundary_edge_dofs_out via ListToMarker().
+       @param[out] dof_edges Optional; local edge index of each DOF.
+       @param[out] dof_boundary_elements Optional; a boundary element containing
+                   each DOF.
+       @param[out] ess_edge_list Optional array of edge indices, in one-to-one
+                                 correspondence with @a ess_tdof_list. An entry
+                                 is -1 when the true DOF is owned by this rank
+                                 but no local edge can be associated with it,
+                                 which can happen for a shared vertex DOF whose
+                                 boundary elements are all on other ranks. */
+   void GetBoundaryLoopEdgeDofs(const Array<int> &boundary_element_indices,
+                                Array<int> &ess_tdof_list,
+                                Array<int> &boundary_edge_dofs_out,
+                                Array<int> *ldof_marker = nullptr,
+                                Array<int> *dof_edges = nullptr,
+                                Array<int> *dof_boundary_elements = nullptr,
+                                Array<int> *ess_edge_list = nullptr);
+
   /** If the given ldof is owned by the current processor, return its local
       tdof number, otherwise return -1 */
   int GetLocalTDofNumber(int ldof) const;
@@ -17,7 +17,6 @@ list(APPEND SRCS
  error.cpp
  gecko.cpp
  globals.cpp
-  glvis_stream.cpp
  hash.cpp
  hash_util.cpp
  isockstream.cpp
@@ -46,7 +45,6 @@ list(APPEND HDRS
  error.hpp
  gecko.hpp
  globals.hpp
-  glvis_stream.hpp
  zstr.hpp
  hash.hpp
  hash_util.hpp
@@ -1108,6 +1108,126 @@ void GroupCommunicator::ReduceEnd(T *ldata, int layout,
   num_requests = 0;
 }

+template <class T>
+void GroupCommunicator::ReduceMarked(T *ldata, const Array<int> &marker,
+                                     int layout,
+                                     void (*Op)(OpData<T>)) const
+{
+   if (comm_lock == 0) { return; }
+   // The above also handles the case (group_buf_size == 0).
+   MFEM_VERIFY(comm_lock == 2, "object is NOT locked for Reduce");
+
+   switch (mode)
+   {
+      case byGroup: // ***** Communication by groups *****
+      {
+         OpData<T> opd;
+         opd.ldata = ldata;
+         Array<int> group_num_req(group_ldof.Size());
+         for (int gr = 1; gr < group_ldof.Size(); gr++)
+         {
+            group_num_req[gr] =
+               gtopo.IAmMaster(gr) ? gtopo.GetGroupSize(gr)-1 : 0;
+         }
+         int idx;
+         while (MPI_Waitany(num_requests, requests, &idx, MPI_STATUS_IGNORE),
+                idx != MPI_UNDEFINED)
+         {
+            int gr = request_marker[idx];
+            if (gr == -1) { continue; } // skip send requests
+
+            // Delay the processing of a group until all receive requests, for
+            // that group, are done:
+            if ((--group_num_req[gr]) != 0) { continue; }
+
+            opd.nldofs = group_ldof.RowSize(gr);
+            // groups without dofs are skipped, so here nldofs > 0.
+
+            opd.buf = (T *)group_buf.GetData() + buf_offsets[gr];
+            opd.ldofs = (layout == 0) ?
+                        group_ldof.GetRow(gr) : group_ltdof.GetRow(gr);
+            opd.nb = gtopo.GetGroupSize(gr)-1;
+
+            // Apply operation only to marked DOFs. The receive buffer is
+            // neighbor-major with stride opd.nldofs, i.e. the contributions to
+            // DOF i are buf[j*opd.nldofs + i] for j = 0 ... opd.nb-1. Setting
+            // nldofs = 1 for a single DOF changes that stride to 1, so the
+            // strided values must first be gathered into a contiguous buffer.
+            Array<T> single_buf(opd.nb);
+            for (int i = 0; i < opd.nldofs; i++)
+            {
+               if (marker[opd.ldofs[i]])
+               {
+                  for (int j = 0; j < opd.nb; j++)
+                  {
+                     single_buf[j] = opd.buf[j*opd.nldofs + i];
+                  }
+
+                  // Create a temporary OpData with just this one DOF
+                  OpData<T> single_opd;
+                  single_opd.ldata = ldata;
+                  single_opd.buf = single_buf.GetData();
+                  single_opd.ldofs = opd.ldofs + i;
+                  single_opd.nldofs = 1;
+                  single_opd.nb = opd.nb;
+
+                  // Apply the operation
+                  Op(single_opd);
+               }
+            }
+         }
+         break;
+      }
+
+      case byNeighbor: // ***** Communication by neighbors *****
+      {
+         MPI_Waitall(num_requests, requests, MPI_STATUSES_IGNORE);
+
+         for (int nbr = 1; nbr < nbr_send_groups.Size(); nbr++)
+         {
+            // In Reduce operation: send_groups <--> recv_groups
+            const int num_recv_groups = nbr_send_groups.RowSize(nbr);
+            if (num_recv_groups > 0)
+            {
+               const int *grp_list = nbr_send_groups.GetRow(nbr);
+               const T *buf = (T*)group_buf.GetData() + buf_offsets[nbr];
+               for (int i = 0; i < num_recv_groups; i++)
+               {
+                  // Custom version of ReduceGroupFromBuffer that checks marker
+                  int gr = grp_list[i];
+                  const int *ldofs = (layout == 0) ?
+                                     group_ldof.GetRow(gr) : group_ltdof.GetRow(gr);
+                  const int nldofs = group_ldof.RowSize(gr);
+
+                  for (int j = 0; j < nldofs; j++)
+                  {
+                     if (marker[ldofs[j]])
+                     {
+                        // Create a temporary OpData with just this one DOF
+                        OpData<T> opd;
+                        opd.ldata = ldata;
+                        opd.buf = const_cast<T*>(buf) + j;
+                        opd.ldofs = ldofs + j;
+                        opd.nldofs = 1;
+                        opd.nb = 1;
+
+                        // Apply the operation
+                        Op(opd);
+                     }
+                  }
+
+                  buf += nldofs;
+               }
+            }
+         }
+         break;
+      }
+   }
+
+   comm_lock = 0; // 0 - no lock
+   num_requests = 0;
+}
+
 template <class T>
 void GroupCommunicator::Sum(OpData<T> opd)
 {
@@ -1171,6 +1291,8 @@ void GroupCommunicator::Max(OpData<T> opd)
 template <class T>
 void GroupCommunicator::BitOR(OpData<T> opd)
 {
+   static_assert(std::is_integral<T>::value,
+                 "BitOR reduction requires an integral type.");
   for (int i = 0; i < opd.nldofs; i++)
   {
      T data = opd.ldata[opd.ldofs[i]];
@@ -1182,6 +1304,33 @@ void GroupCommunicator::BitOR(OpData<T> opd)
   }
 }

+template <class T>
+void GroupCommunicator::MaxAbs(OpData<T> opd)
+{
+   for (int i = 0; i < opd.nldofs; i++)
+   {
+      T data = opd.ldata[opd.ldofs[i]];
+      T abs_data = std::abs(data);
+
+      for (int j = 0; j < opd.nb; j++)
+      {
+         T b = opd.buf[j*opd.nldofs+i];
+         T abs_b = std::abs(b);
+
+         // On an equal-magnitude tie keep the more positive value, so
+         // opposite-sign ties resolve deterministically to the positive one.
+         if (abs_data < abs_b || (abs_data == abs_b && data < b))
+         {
+            data = b;
+            abs_data = abs_b;
+         }
+      }
+
+      opd.ldata[opd.ldofs[i]] = data;
+   }
+}
+
+
 void GroupCommunicator::PrintInfo(std::ostream &os) const
 {
   char c = '\0';
@@ -1318,18 +1467,24 @@ template void GroupCommunicator::BcastEnd<int>(int *, int) const;
 template void GroupCommunicator::ReduceBegin<int>(const int *) const;
 template void GroupCommunicator::ReduceEnd<int>(
   int *, int, void (*)(OpData<int>)) const;
+template void GroupCommunicator::ReduceMarked<int>(
+   int*, const Array<int>&, int, void (*)(OpData<int>)) const;

 template void GroupCommunicator::BcastBegin<double>(double *, int) const;
 template void GroupCommunicator::BcastEnd<double>(double *, int) const;
 template void GroupCommunicator::ReduceBegin<double>(const double *) const;
 template void GroupCommunicator::ReduceEnd<double>(
   double *, int, void (*)(OpData<double>)) const;
+template void GroupCommunicator::ReduceMarked<double>(
+   double*, const Array<int>&, int, void (*)(OpData<double>)) const;

 template void GroupCommunicator::BcastBegin<float>(float *, int) const;
 template void GroupCommunicator::BcastEnd<float>(float *, int) const;
 template void GroupCommunicator::ReduceBegin<float>(const float *) const;
 template void GroupCommunicator::ReduceEnd<float>(
   float *, int, void (*)(OpData<float>)) const;
+template void GroupCommunicator::ReduceMarked<float>(
+   float*, const Array<int>&, int, void (*)(OpData<float>)) const;

 // @endcond

@@ -1338,14 +1493,17 @@ template void GroupCommunicator::Sum<int>(OpData<int>);
 template void GroupCommunicator::Min<int>(OpData<int>);
 template void GroupCommunicator::Max<int>(OpData<int>);
 template void GroupCommunicator::BitOR<int>(OpData<int>);
+template void GroupCommunicator::MaxAbs<int>(OpData<int>);

 template void GroupCommunicator::Sum<double>(OpData<double>);
 template void GroupCommunicator::Min<double>(OpData<double>);
 template void GroupCommunicator::Max<double>(OpData<double>);
+template void GroupCommunicator::MaxAbs<double>(OpData<double>);

 template void GroupCommunicator::Sum<float>(OpData<float>);
 template void GroupCommunicator::Min<float>(OpData<float>);
 template void GroupCommunicator::Max<float>(OpData<float>);
+template void GroupCommunicator::MaxAbs<float>(OpData<float>);


 #ifdef __bgq__
@@ -22,6 +22,7 @@
 #include "globals.hpp"
 #include <mpi.h>
 #include <cstdint>
+#include <type_traits>

 // can't directly use MPI_CXX_BOOL because Microsoft's MPI implementation
 // doesn't include MPI_CXX_BOOL. Fallback to MPI_C_BOOL if unavailable.
@@ -408,14 +409,38 @@ public:
   template <class T> void Reduce(Array<T> &ldata, void (*Op)(OpData<T>)) const
   { Reduce<T>((T *)ldata, Op); }

-   /// Reduce operation Sum, instantiated for int and double
+   /// Reduce operation Sum, instantiated for int, double and float
   template <class T> static void Sum(OpData<T>);
-   /// Reduce operation Min, instantiated for int and double
+   /// Reduce operation Min, instantiated for int, double and float
   template <class T> static void Min(OpData<T>);
-   /// Reduce operation Max, instantiated for int and double
+   /// Reduce operation Max, instantiated for int, double and float
   template <class T> static void Max(OpData<T>);
   /// Reduce operation bitwise OR, instantiated for int only
   template <class T> static void BitOR(OpData<T>);
+   /// Reduce operation selecting the signed value with the largest absolute
+   /// value, instantiated for int, double and float. The result keeps its sign;
+   /// it is not the non-negative absolute value. Equal-magnitude ties are
+   /// broken deterministically toward the more positive value, so opposite-sign
+   /// ties resolve to the positive one regardless of accumulation order.
+   template <class T> static void MaxAbs(OpData<T>);
+
+   /** @brief Finalize reduction operation started with ReduceBegin(), but only apply
+       the reduction to DOFs marked in the marker array.
+
+       @note The reduction is carried out in the signed type @a T, so the result
+       is signed even for bitwise operations.
+   */
+   template <class T>
+   void ReduceMarked(T *ldata, const Array<int> &marker, int layout,
+                     void (*Op)(OpData<T>)) const;
+
+   /** @brief Reduce within each group where the master is the root, but only for marked DOFs. */
+   template <class T>
+   void Reduce(T *ldata, const Array<int> &marker, void (*Op)(OpData<T>)) const
+   {
+      ReduceBegin(ldata);
+      ReduceMarked(ldata, marker, 0, Op);
+   }

   /// Print information about the GroupCommunicator from all MPI ranks.
   void PrintInfo(std::ostream &out = mfem::out) const;
@@ -1,188 +0,0 @@
-// 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.
-
-#include <cassert>
-#include <istream>
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_GLVIS
-
-#include "glvis_stream.hpp"
-
-#ifdef MFEM_USE_MPI
-#include <mpi.h>
-#include <limits>
-#include <numeric>
-#endif
-
-#include "../fem/geom.hpp"
-thread_local mfem::GeometryRefiner GLVisGeometryRefiner;
-
-// Use local declaration to avoid circular dependency when fetching GLVis
-extern int GLVisStreamSession(
-   bool fix_elem_orient,
-   bool save_coloring,
-   bool keep_attr,
-   bool headless,
-   const std::string& plot_caption,
-   const std::string& data_type,
-   std::vector<std::unique_ptr<std::istream>>&& streams);
-
-namespace mfem
-{
-
-#ifdef MFEM_USE_MPI
-namespace
-{
-glvis_data MakeGlVisData()
-{
-   int size, rank;
-   MPI_Comm_size(MPI_COMM_WORLD, &size);
-   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-   return glvis_data(size == 1, size, rank == 0);
-}
-} // namespace
-#endif
-
-glvis_stream::glvis_stream(): std::iostream(nullptr),
-#ifdef MFEM_USE_MPI
-   data(MakeGlVisData())
-#else
-   data(true, 1, true)
-#endif
-{
-   std::iostream::rdbuf(data.stream.rdbuf());
-}
-
-glvis_stream& glvis_stream::operator<<(ostream_manipulator pf)
-{
-   pf(static_cast<std::ostream&>(*this));
-   this->flush();
-   this->operator()();
-   return *this;
-}
-
-void glvis_stream::operator()()
-{
-   if (data.serial)
-   {
-      const auto size = this->size();
-      data.offsets.resize(2);
-      data.offsets[0] = 0, data.offsets[1] = size;
-      data.total_size = size;
-   }
-   else
-   {
-      serialize();
-   }
-
-   this->reset(); // reset the local buffer for reuse
-
-   if (data.mpi_root)
-   {
-      MFEM_VERIFY(data.mpi_size >= 0 &&
-                  (size_t) data.mpi_size == data.offsets.size() - 1,
-                  "Invalid MPI size");
-
-      data.streams.clear();
-      data.type.clear();
-
-      // loop over all input streams
-      for (int k = 0; k < data.mpi_size; ++k)
-      {
-         const size_t offset = data.offsets[k];
-         const size_t size = data.offsets[k+1] - data.offsets[k];
-
-         // add a new stream for this rank's data
-         data.streams.emplace_back(std::make_unique<std::stringstream>());
-         data.streams.back()->write(data.stream.str().data() + offset, size);
-
-         auto stream = data.streams.back().get();
-         if (!(*stream)) { break; }
-
-         *stream >> std::ws >> data.type >> std::ws;
-
-         if (data.type == "parallel") // Handle parallel data
-         {
-            int is_mpi_size, is_mpi_rank;
-            *stream >> is_mpi_size >> is_mpi_rank;
-            assert(is_mpi_size == static_cast<int>(data.mpi_size));
-            assert(is_mpi_rank == static_cast<int>(k));
-         }
-         else if (data.type != "mesh" && data.type != "solution")
-         {
-            MFEM_ABORT("Stream: unknown command: " << data.type);
-         }
-      }
-   }
-
-   if (!data.mpi_root) { return; }
-
-   constexpr bool fix_elem_orien = true;
-   constexpr bool save_coloring = true;
-   constexpr bool keep_attr = false;
-   constexpr bool headless = false;
-   const std::string plot_caption {};
-   std::vector<std::unique_ptr<std::istream>> istreams;
-   istreams.reserve(data.streams.size());
-   for (auto &s : data.streams) { istreams.push_back(std::move(s)); }
-   GLVisStreamSession(fix_elem_orien,
-                      save_coloring,
-                      keep_attr,
-                      headless,
-                      plot_caption,
-                      data.type,
-                      std::move(istreams));
-}
-
-void glvis_stream::serialize()
-{
-#ifdef MFEM_USE_MPI
-   const std::string local = data.stream.str();
-   MFEM_VERIFY(local.size() <= static_cast<size_t>
-               (std::numeric_limits<int>::max()),
-               "GLVis stream is too large for MPI_Gatherv");
-   const int local_size = static_cast<int>(local.size());
-
-   std::vector<int> sizes(data.mpi_size);
-   MFEM_VERIFY(MPI_Allgather(&local_size, 1, MPI_INT,
-                             sizes.data(), 1, MPI_INT,
-                             MPI_COMM_WORLD) == MPI_SUCCESS,
-               "MPI_Allgather failed");
-
-   if (data.mpi_root)
-   {
-      data.offsets.resize(data.mpi_size + 1);
-      data.offsets[0] = 0;
-      std::partial_sum(sizes.begin(), sizes.end(), data.offsets.begin() + 1);
-      data.total_size = data.offsets[data.mpi_size];
-   }
-
-   std::vector<char> recvbuf(data.mpi_root ? data.total_size : 0);
-   MFEM_VERIFY(MPI_Gatherv(local.data(), local_size, MPI_CHAR,
-                           data.mpi_root ? recvbuf.data() : nullptr,
-                           data.mpi_root ? sizes.data() : nullptr,
-                           data.mpi_root ? data.offsets.data() : nullptr,
-                           MPI_CHAR, 0, MPI_COMM_WORLD) == MPI_SUCCESS,
-               "MPI_Gatherv failed");
-
-   if (data.mpi_root)
-   {
-      reset();
-      data.stream.write(recvbuf.data(), data.total_size);
-   }
-#endif // MFEM_USE_MPI
-}
-
-} // namespace mfem
-
-#endif // MFEM_USE_GLVIS
@@ -1,83 +0,0 @@
-// 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.
-#pragma once
-
-#include <memory>
-#include <sstream>
-#include <string>
-#include <vector>
-
-namespace mfem
-{
-
-struct glvis_data
-{
-   const bool serial;
-   const int mpi_size;
-   const bool mpi_root;
-   std::stringstream stream;
-   std::vector<std::unique_ptr<std::stringstream>> streams;
-   int total_size;
-   std::vector<int> offsets;
-   std::string type;
-
-   glvis_data(const bool serial, const int size, const bool root):
-      serial(serial), mpi_size(size), mpi_root(root),
-      total_size(0), type({}) {}
-};
-
-class glvis_stream: public std::iostream
-{
-   glvis_data data;
-   void serialize();
-
-public:
-   glvis_stream();
-
-   glvis_stream(glvis_stream &&) = delete;
-   glvis_stream(const glvis_stream &) = delete;
-   glvis_stream &operator=(const glvis_stream &) = delete;
-   glvis_stream &operator=(glvis_stream &&) = delete;
-
-   ~glvis_stream() = default;
-
-   size_t size() { return data.stream.tellp(); }
-   std::streamsize precision() const { return std::iostream::precision(); }
-   std::streamsize precision(std::streamsize new_prec)
-   { return std::iostream::precision(new_prec); }
-
-   using ostream_manipulator = std::ostream& (*)(std::ostream&);
-   glvis_stream& operator<<(ostream_manipulator pf);
-
-   template<typename T>
-   glvis_stream& operator<<(const T& val)
-   {
-      static_cast<std::ostream&>(*this) << val;
-      return *this;
-   }
-
-   int open(const char *, int) { return 0; }
-   bool is_open() const { return true; }
-   int close() { return 0; }
-
-   void flush() { std::iostream::flush(); }
-
-   void reset()
-   {
-      data.stream.clear();
-      data.stream.seekg(0, std::ios::beg);
-      data.stream.seekp(0, std::ios::beg);
-   }
-
-   void operator()();
-};
-
-} // namespace mfem
@@ -1087,4 +1087,3 @@ socketstream::~socketstream()
 }

 } // namespace mfem
-
@@ -13,7 +13,6 @@
 #define MFEM_SOCKETSTREAM

 #include "../config/config.hpp"
-
 #include "error.hpp"
 #include "globals.hpp"

@@ -100,9 +100,6 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_GSLIB
      "MFEM_USE_GSLIB\n"
 #endif
-#ifdef MFEM_USE_GLVIS
-      "MFEM_USE_GLVIS\n"
-#endif
 #ifdef MFEM_USE_HDF5
      "MFEM_USE_HDF5\n"
 #endif
@@ -119,8 +119,7 @@ $(if $(word 2,$(SRC)),$(error Spaces in SRC = "$(SRC)" are not supported))
 MFEM_GIT_STRING = $(shell [ -d $(MFEM_DIR)/.git ] && git -C $(MFEM_DIR) \
   describe --all --long --abbrev=40 --dirty --always 2> /dev/null)

-EXAMPLE_SUBDIRS = amgx caliper ginkgo glvis hiop petsc pumi sundials \
- superlu moonolith
+EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples

@@ -164,8 +163,6 @@ MFEM_BUILD_DIR := $(BUILD_DIR)

 CONFIG_MK = $(BLD)config/config.mk

-GLVIS_MK = $(SRC)config/glvis.mk
-
 DEFAULTS_MK = $(SRC)config/defaults.mk
 include $(DEFAULTS_MK)

@@ -294,11 +291,6 @@ ifeq ($(MFEM_USE_HIP),YES)
   endif
 endif

-# GLVis configuration
-ifeq ($(MFEM_USE_GLVIS),YES)
-	GLVIS_DIR:=$(abspath $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(GLVIS_DIR)))
-endif
-
 DEP_CXX ?= $(MFEM_CXX)

 # Check legacy OpenMP configuration
@@ -315,7 +307,7 @@ endif
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
 SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
 GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
- TRIBOL BENCHMARK MOONOLITH ALGOIM CUDSS GLVIS
+ TRIBOL BENCHMARK MOONOLITH ALGOIM CUDSS


 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
@@ -385,8 +377,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
 MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
 MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL MFEM_USE_ALGOIM MFEM_USE_ENZYME\
 MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE\
- MFEM_USE_CUDSS MFEM_CUDSS_COMM_LIB MFEM_CUDSS_THREADING_LIB\
- MFEM_USE_GLVIS
+ MFEM_USE_CUDSS MFEM_CUDSS_COMM_LIB MFEM_CUDSS_THREADING_LIB

 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -486,8 +477,7 @@ OKL_DIRS = fem
 %:	%.cpp

 # Default rule.
-lib: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT)) \
-	$(if $(filter YES,$(MFEM_USE_GLVIS)),$(if $(static), $(GLVIS_DIR)/lib/libglvis.a))
+lib: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT))

 # Flags used for compiling all source files.
 MFEM_BUILD_FLAGS = $(MFEM_PICFLAG) $(MFEM_CPPFLAGS) $(MFEM_CXXFLAGS)\
@@ -520,14 +510,6 @@ $(BLD)libmfem.$(SO_EXT): $(BLD)libmfem.$(SO_VER)
 	cd $(@D) && ln -sf $(<F) $(@F)
 	@$(MAKE) deprecation-warnings

-ifeq ($(MFEM_USE_GLVIS),YES)
-$(GLVIS_DIR)/lib/libglvis.a: $(BLD)libmfem.a
-	$(if $(wildcard $(GLVIS_DIR)/makefile),,$(error No makefile in GLVIS_DIR: $(GLVIS_DIR)))
-	@$(MAKE) -C $(GLVIS_DIR) -j $(shell getconf _NPROCESSORS_ONLN 2>/dev/null || echo 1) \
-		MFEM_DIR=$(BUILD_REAL_DIR) \
-		GLVIS_USE_LOGO=NO GLVIS_USE_LIBPNG=YES lib/libglvis.a
-endif
-
 # If some of the external libraries are build without -fPIC, linking shared MFEM
 # library may fail. In such cases, one may set EXT_LIBS on the command line.
 EXT_LIBS = $(MFEM_EXT_LIBS)
@@ -622,7 +604,6 @@ clean: $(addsuffix /clean,$(EM_DIRS) $(TEST_DIRS))

 distclean: clean config/clean doc/clean
 	rm -rf mfem/
-	$(if $(filter YES,$(MFEM_USE_GLVIS)),-$(MAKE) -C $(GLVIS_DIR) distclean)

 # User-definable install permissions.
 # Install permissions for everything except directories and binaries:
@@ -647,14 +628,10 @@ INSTALL_SHARED_LIB = $(MFEM_CXX) $(MFEM_LINK_FLAGS) $(INSTALL_SOFLAGS)\
   cd $(PREFIX_LIB) && chmod $(INSTALL_BIN_PERM) libmfem.$(SO_VER) && \
   ( umask $(INSTALLMASK) && ln -sf libmfem.$(SO_VER) libmfem.$(SO_EXT) )

-install: $(if $(static),$(BLD)libmfem.a) \
-			$(if $(shared),$(BLD)libmfem.$(SO_EXT)) \
-			$(if $(filter YES,$(MFEM_USE_GLVIS)),\
-				$(if $(static),$(GLVIS_DIR)/lib/libglvis.a))
+install: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT))
 	$(MKINSTALLDIR) $(PREFIX_LIB)
 # install static and/or shared library
 	$(if $(static),$(INSTALLDEF) $(BLD)libmfem.a $(PREFIX_LIB))
-	$(if $(filter YES,$(MFEM_USE_GLVIS)),$(if $(static),$(INSTALLDEF) $(GLVIS_DIR)/lib/libglvis.a) $(PREFIX_LIB))
 	$(if $(shared),$(INSTALL_SHARED_LIB))
 # install top level includes
 	$(MKINSTALLDIR) $(PREFIX_INC)/mfem
@@ -801,7 +778,6 @@ status info:
 	$(info MFEM_USE_PARELAG       = $(MFEM_USE_PARELAG))
 	$(info MFEM_USE_TRIBOL        = $(MFEM_USE_TRIBOL))
 	$(info MFEM_USE_ENZYME        = $(MFEM_USE_ENZYME))
-	$(info MFEM_USE_GLVIS         = $(MFEM_USE_GLVIS))
 	$(info MFEM_CXX               = $(value MFEM_CXX))
 	$(info MFEM_HOST_CXX          = $(value MFEM_HOST_CXX))
 	$(info MFEM_CPPFLAGS          = $(value MFEM_CPPFLAGS))
@@ -30,9 +30,6 @@
 #ifdef MFEM_USE_ADIOS2
 #include "general/adios2stream.hpp"
 #endif // MFEM_USE_ADIOS2
-#ifdef MFEM_USE_GLVIS
-#include "general/glvis_stream.hpp"
-#endif // MFEM_USE_GLVIS
 #include "general/isockstream.hpp"
 #include "general/osockstream.hpp"
 #include "general/socketstream.hpp"
@@ -97,6 +97,7 @@ set(UNIT_TESTS_SRCS
  fem/test_2d_bilininteg.cpp
  fem/test_3d_bilininteg.cpp
  fem/test_assembly_levels.cpp
+  fem/test_bdr_edgedof.cpp
  fem/test_bilinearform.cpp
  fem/test_block_operators.cpp
  fem/test_blocknonlinearform.cpp
@@ -0,0 +1,705 @@
+// 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.
+
+#include "unit_tests.hpp"
+#include "mfem.hpp"
+#include "../mesh/mesh_test_utils.hpp"
+
+#include <set>
+#include <unordered_set>
+#include <vector>
+
+using namespace mfem;
+
+#ifdef MFEM_USE_MPI
+
+TEST_CASE("BoundaryEdgeDOFsPartitionInvariant",
+          "[Parallel][ParMesh][BoundaryEdgeDOFs]")
+{
+   constexpr int orientation = 3;
+   constexpr int order = 1;
+
+   // Use all available MPI processes for partitioning
+   const int test_num_procs = Mpi::WorldSize();
+
+   // Create base mesh
+   Mesh base_mesh = OrientedTriFaceMesh(orientation, true);
+   base_mesh.UniformRefinement();
+   const int n_elements = base_mesh.GetNE();
+
+   // Use a small set of representative partitionings
+   std::vector<std::vector<int>> all_partitionings;
+   // 1. All elements on rank 0
+   all_partitionings.push_back(std::vector<int>(n_elements, 0));
+   if (test_num_procs > 1)
+   {
+      // 2. Block partition: first half on rank 0, second half on last rank
+      std::vector<int> &block = all_partitionings.emplace_back(n_elements);
+      for (int i = 0; i < n_elements; i++)
+      {
+         block[i] = (i < n_elements/2) ? 0 : test_num_procs-1;
+      }
+
+      // 3. Round-robin partition: elements assigned cyclically to all ranks
+      std::vector<int> &round_robin = all_partitionings.emplace_back(n_elements);
+      for (int i = 0; i < n_elements; i++)
+      {
+         round_robin[i] = i % test_num_procs;
+      }
+   }
+
+   // Create reusable FEC
+   ND_FECollection fec(order, 3);
+
+   std::vector<int> all_results;
+   all_results.reserve(all_partitionings.size());
+
+   // Test each partitioning
+   for (const auto& partition : all_partitionings)
+   {
+      // Create parallel mesh with current partitioning
+      Mesh test_mesh = OrientedTriFaceMesh(orientation, true);
+      test_mesh.UniformRefinement();
+      // For single process, use default partitioning; for multiple, use custom partition
+      ParMesh pmesh = (test_num_procs == 1) ?
+                      ParMesh(MPI_COMM_WORLD, test_mesh) :
+                      ParMesh(MPI_COMM_WORLD, test_mesh, partition.data());
+
+      // Create finite element space
+      ParFiniteElementSpace fespace(&pmesh, &fec);
+
+      // Extract boundary edge DOFs
+      Array<int> ess_tdof_list;
+      Array<int> boundary_edge_ldofs;
+      std::vector<Array<int>> attr_to_elements;
+
+      // Select the shared face to be the tested boundary
+      int bdr_attr = pmesh.bdr_attributes.Max();
+      Array<int> bdr_attrs(1);
+      bdr_attrs[0] = bdr_attr;
+
+      fespace.GetBoundaryElementsByAttribute(bdr_attrs, attr_to_elements);
+      Array<int> boundary_elements = attr_to_elements[0];
+
+      Array<int> dof_edges, dof_boundary_elements, ess_edge_list;
+
+      fespace.GetBoundaryLoopEdgeDofs(boundary_elements, ess_tdof_list,
+                                      boundary_edge_ldofs, nullptr, &dof_edges,
+                                      &dof_boundary_elements, &ess_edge_list);
+
+      // Collect total boundary edge DOFs
+      int local_dofs = boundary_edge_ldofs.Size();
+      int total_dofs;
+      MPI_Allreduce(&local_dofs, &total_dofs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+
+      all_results.push_back(total_dofs);
+   }
+
+   // The set of boundary edge DOFs is a property of the mesh geometry and must
+   // not depend on how the elements are distributed across ranks. Each result
+   // is the global count of selected boundary edge DOFs for one partitioning, so
+   // if the method correctly removes the artificial edges introduced at
+   // processor boundaries, every partitioning yields the same total. A mismatch
+   // means some partition kept or dropped a DOF that another did not.
+   REQUIRE(!all_results.empty());
+   // One refinement splits the triangular face into four sub-triangles. Its
+   // perimeter has six loop edges (order-1 ND: one DOF per edge); the three
+   // interior edges of the middle sub-triangle are shared and correctly dropped.
+   constexpr int expected = 6;
+   for (int result : all_results)
+   {
+      REQUIRE(result == expected);
+   }
+}
+
+TEST_CASE("BoundaryEdgeDOFsBasicFunctionality",
+          "[Parallel][ParMesh][BoundaryEdgeDOFs]")
+{
+   const int orientation = GENERATE(1, 3, 5);
+   const int order = GENERATE(1, 2);
+
+   CAPTURE(orientation, order);
+
+   // Create test mesh
+   Mesh mesh = OrientedTriFaceMesh(orientation, true);
+   mesh.UniformRefinement();
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+
+   // Create finite element space
+   ND_FECollection fec(order, 3);
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+
+   // Test boundary edge DOF extraction
+   Array<int> ess_tdof_list;
+   Array<int> boundary_edge_ldofs;
+   Array<int> ldof_marker;
+   std::vector<Array<int>> attr_to_elements;
+
+   // Get boundary elements for the shared face
+   int bdr_attr = pmesh.bdr_attributes.Max();
+   Array<int> bdr_attrs(1);
+   bdr_attrs[0] = bdr_attr;
+
+   fespace.GetBoundaryElementsByAttribute(bdr_attrs, attr_to_elements);
+   Array<int> boundary_elements = attr_to_elements[0];
+
+   Array<int> dof_edges, dof_boundary_elements, ess_edge_list;
+
+   fespace.GetBoundaryLoopEdgeDofs(boundary_elements, ess_tdof_list,
+                                   boundary_edge_ldofs, &ldof_marker, &dof_edges,
+                                   &dof_boundary_elements, &ess_edge_list);
+
+   // Basic validation
+   REQUIRE(ldof_marker.Size() == fespace.GetVSize());
+   REQUIRE(ess_tdof_list.Size() >= 0);
+
+   // The output arrays share a single indexing, so they must have equal size.
+   REQUIRE(boundary_edge_ldofs.Size() == dof_edges.Size());
+   REQUIRE(dof_edges.Size() == dof_boundary_elements.Size());
+
+   // Verify all boundary edge DOFs are marked in ldof_marker
+   for (int dof : boundary_edge_ldofs)
+   {
+      REQUIRE(ldof_marker[dof] == 1);
+   }
+}
+
+// Helper function to compute boundary loop length
+real_t ComputeBoundaryLoopLength(ParMesh* pmesh, const Array<int>& dof_edges)
+{
+   real_t local_length = 0.0;
+   std::unordered_set<int> processed_edges;
+
+   for (int i = 0; i < dof_edges.Size(); i++)
+   {
+      int edge_id = dof_edges[i];
+      if (!processed_edges.insert(edge_id).second) { continue; }
+
+      Array<int> edge_verts;
+      pmesh->GetEdgeVertices(edge_id, edge_verts);
+      const real_t* v0 = pmesh->GetVertex(edge_verts[0]);
+      const real_t* v1 = pmesh->GetVertex(edge_verts[1]);
+
+      real_t edge_length = 0.0;
+      for (int d = 0; d < pmesh->SpaceDimension(); d++)
+      {
+         real_t diff = v1[d] - v0[d];
+         edge_length += diff * diff;
+      }
+      local_length += sqrt(edge_length);
+   }
+   return local_length;
+}
+
+TEST_CASE("BoundaryEdgeDOFsNestedCubes",
+          "[Parallel][ParMesh][BoundaryEdgeDOFs]")
+{
+   const int order = GENERATE(1, 2);
+
+   // Expected processor-invariant results for nested cubes mesh (1 refinement)
+   // order=1: 16 tdofs, sum=16.0, length=2.0
+   // order=2: 32 tdofs, sum=32.0, length=2.0
+   int exp_tdofs = (order == 1) ? 16 : 32;
+   real_t exp_sum = (order == 1) ? real_t(16.0) : real_t(32.0);
+   real_t exp_length = real_t(2.0);
+
+   struct BoundaryTest
+   {
+      int attr_value;
+      Vector normal;
+      std::string name;
+   };
+
+   std::vector<BoundaryTest> boundary_tests =
+   {
+      {7, Vector({0, 0, -1}), "-z"},
+      {8, Vector({0, 0, 1}), "+z"},
+      {9, Vector({0, -1, 0}), "-y"},
+      {10, Vector({1, 0, 0}), "+x"},
+      {11, Vector({0, 1, 0}), "+y"},
+      {12, Vector({-1, 0, 0}), "-x"}
+   };
+
+   const char* mesh_file = "../../data/nested_cubes.msh";
+   Mesh mesh(mesh_file, 1, 1);
+   mesh.UniformRefinement();
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+
+   ND_FECollection fec(order, 3);
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+
+   int num_procs;
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+
+   for (const auto& test : boundary_tests)
+   {
+      CAPTURE(test.name, test.attr_value, order, num_procs);
+
+      std::vector<Array<int>> attr_to_elements;
+      Array<int> bdr_attrs(1);
+      bdr_attrs[0] = test.attr_value;
+
+      fespace.GetBoundaryElementsByAttribute(bdr_attrs, attr_to_elements);
+      Array<int> boundary_elements = attr_to_elements[0];
+
+      Array<int> ess_tdof_list;
+      Array<int> ldof_marker;
+      Array<int> boundary_edge_ldofs;
+      Array<int> dof_edges, dof_boundary_elements, ess_edge_list;
+
+      fespace.GetBoundaryLoopEdgeDofs(boundary_elements, ess_tdof_list,
+                                      boundary_edge_ldofs, &ldof_marker, &dof_edges,
+                                      &dof_boundary_elements, &ess_edge_list);
+
+      Array<int> dof_orientations;
+      fespace.ComputeLoopEdgeOrientations(dof_edges, dof_boundary_elements,
+                                          test.normal, dof_orientations);
+
+      ParGridFunction x(&fespace);
+      x = real_t(0.0);
+      for (int i = 0; i < boundary_edge_ldofs.Size(); i++)
+      {
+         x(boundary_edge_ldofs[i]) = real_t(1.0) * dof_orientations[i];
+      }
+
+      GroupCommunicator *gc = fespace.ScalarGroupComm();
+      Array<int> global_marker(ldof_marker);
+      gc->Reduce<int>(global_marker.GetData(), GroupCommunicator::BitOR<int>);
+      gc->Bcast(global_marker);
+
+      Array<real_t> values(x.GetData(), x.Size());
+      gc->ReduceBegin(values.GetData());
+      gc->ReduceMarked<real_t>(values.GetData(), global_marker, 0,
+                               GroupCommunicator::MaxAbs<real_t>);
+      gc->Bcast(values.GetData());
+      delete gc;
+
+      Vector x_true;
+      x.GetTrueDofs(x_true);
+
+      int local_nonzero_tdofs = 0;
+      real_t local_tdof_sum = 0.0;
+      for (int tdof = 0; tdof < x_true.Size(); tdof++)
+      {
+         real_t tdof_value = x_true(tdof);
+         if (abs(tdof_value) > 1e-12)
+         {
+            local_nonzero_tdofs++;
+            local_tdof_sum += abs(tdof_value);
+         }
+      }
+
+      real_t local_length = ComputeBoundaryLoopLength(&pmesh, dof_edges);
+
+      int global_nonzero_tdofs;
+      real_t global_tdof_sum, total_length;
+      MPI_Allreduce(&local_nonzero_tdofs, &global_nonzero_tdofs, 1, MPI_INT, MPI_SUM,
+                    MPI_COMM_WORLD);
+      MPI_Allreduce(&local_tdof_sum, &global_tdof_sum, 1,
+                    MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                    MPI_COMM_WORLD);
+      MPI_Allreduce(&local_length, &total_length, 1,
+                    MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                    MPI_COMM_WORLD);
+
+      // Verify processor-invariant results match expected values
+      REQUIRE(global_nonzero_tdofs == exp_tdofs);
+      REQUIRE(abs(global_tdof_sum - exp_sum) < real_t(1e-12));
+      REQUIRE(abs(total_length - exp_length) < real_t(1e-12));
+   }
+}
+
+TEST_CASE("BoundaryEdgeDOFs2DSquareInSquare",
+          "[Parallel][ParMesh][BoundaryEdgeDOFs]")
+{
+   // Test 2D boundary edge DOF extraction using square-in-square mesh
+   constexpr int order = 2;
+
+   // Test multiple inner boundary attributes
+   std::vector<int> inner_attrs_to_test = {5, 6, 7, 8};
+
+   // Load 2D square-in-square mesh from file
+   const char* mesh_file = "../../data/square_in_square.msh";
+   Mesh serial_mesh(mesh_file, 1, 1);
+   serial_mesh.UniformRefinement();
+
+   int num_procs = Mpi::WorldSize();
+
+   // Test each boundary attribute
+   for (int inner_attr : inner_attrs_to_test)
+   {
+      CAPTURE(inner_attr); // Capture the attribute being tested for better test output
+
+      // Test that results are consistent across different mesh partitionings
+      const int n_elements = serial_mesh.GetNE();
+
+      // Generate multiple different partitionings
+      std::vector<std::vector<int>> all_partitionings;
+
+      // 1. All elements on rank 0
+      all_partitionings.push_back(std::vector<int>(n_elements, 0));
+
+      if (num_procs > 1)
+      {
+         // 2. Block partition: first half on rank 0, second half on last rank
+         std::vector<int> block(n_elements);
+         for (int i = 0; i < n_elements; i++)
+         {
+            block[i] = (i < n_elements/2) ? 0 : num_procs-1;
+         }
+         all_partitionings.push_back(block);
+
+         // 3. Round-robin partition: elements assigned cyclically to all ranks
+         std::vector<int> round_robin(n_elements);
+         for (int i = 0; i < n_elements; i++)
+         {
+            round_robin[i] = i % num_procs;
+         }
+         all_partitionings.push_back(round_robin);
+      }
+
+      ND_FECollection fec(order, 2);
+      std::vector<int> all_dof_results;
+      all_dof_results.reserve(all_partitionings.size());
+
+      // Test each partitioning
+      for (const auto& partition : all_partitionings)
+      {
+         // Create parallel mesh with current partitioning
+         Mesh test_mesh(mesh_file, 1, 1);
+         test_mesh.UniformRefinement();
+         ParMesh pmesh = (num_procs == 1) ?
+                         ParMesh(MPI_COMM_WORLD, test_mesh) :
+                         ParMesh(MPI_COMM_WORLD, test_mesh, partition.data());
+
+         ParFiniteElementSpace fespace(&pmesh, &fec);
+
+         // Find boundary elements with the inner attribute
+         std::vector<Array<int>> attr_to_elements;
+         Array<int> inner_attrs(1);
+         inner_attrs[0] = inner_attr;
+         fespace.GetBoundaryElementsByAttribute(inner_attrs, attr_to_elements);
+
+         Array<int> inner_boundary_elements = attr_to_elements[0];
+
+         Array<int> ess_tdofs, ess_edges;
+         Array<int> boundary_dofs;
+         Array<int> dof_edges, dof_boundary_elements;
+
+         fespace.GetBoundaryLoopEdgeDofs(inner_boundary_elements, ess_tdofs,
+                                         boundary_dofs, nullptr, &dof_edges,
+                                         &dof_boundary_elements, &ess_edges);
+
+         // The output arrays share one indexing, so their sizes must match.
+         REQUIRE(boundary_dofs.Size() == dof_edges.Size());
+         REQUIRE(dof_edges.Size() == dof_boundary_elements.Size());
+
+         // Gather global counts for this partitioning
+         int local_dof_count = boundary_dofs.Size();
+         int global_dof_count;
+         MPI_Allreduce(&local_dof_count, &global_dof_count, 1, MPI_INT, MPI_SUM,
+                       MPI_COMM_WORLD);
+
+         all_dof_results.push_back(global_dof_count);
+      }
+
+      // Verify all partitionings give identical results
+      REQUIRE(!all_dof_results.empty());
+
+      int expected_dofs = all_dof_results[0];
+      for (int result : all_dof_results)
+      {
+         REQUIRE(result == expected_dofs);
+      }
+   } // End of inner_attr loop
+}
+
+TEST_CASE("BoundaryEdgeDOFsSharedDOFsAreOwnedBySomeRank",
+          "[Parallel][ParMesh][BoundaryEdgeDOFs]")
+{
+   // Every selected shared DOF must appear in exactly one rank's ess_tdof_list.
+   // Only the group master owns the corresponding true DOF and returns a
+   // non-negative value from GetLocalTDofNumber(), so if the master holds none
+   // of the selected boundary elements the DOF would be emitted by no rank at
+   // all unless the local marker is synchronized across the sharing group.
+   const int nranks = Mpi::WorldSize();
+   if (nranks < 2) { return; }
+
+   constexpr int order = 1;
+   ND_FECollection fec(order, 3);
+
+   for (int orientation : {1, 3, 5})
+   {
+      Mesh probe = OrientedTriFaceMesh(orientation, true);
+      probe.UniformRefinement();
+      const int ne = probe.GetNE();
+
+      // Several partitionings, to vary which rank masters each shared group
+      std::vector<std::vector<int>> partitionings;
+      {
+         std::vector<int> round_robin(ne), block(ne), strided(ne);
+         for (int i = 0; i < ne; i++)
+         {
+            round_robin[i] = i % nranks;
+            block[i] = (i < ne/2) ? 0 : nranks-1;
+            strided[i] = (i * 7 + 3) % nranks;
+         }
+         partitionings = {round_robin, block, strided};
+      }
+
+      for (const auto &partition : partitionings)
+      {
+         Mesh mesh = OrientedTriFaceMesh(orientation, true);
+         mesh.UniformRefinement();
+         ParMesh pmesh(MPI_COMM_WORLD, mesh, partition.data());
+         ParFiniteElementSpace fes(&pmesh, &fec);
+
+         const int bdr_attr = pmesh.bdr_attributes.Max();
+         Array<int> bdr_attrs(1);
+         bdr_attrs[0] = bdr_attr;
+         std::vector<Array<int>> attr_to_elements;
+         fes.GetBoundaryElementsByAttribute(bdr_attrs, attr_to_elements);
+         Array<int> bdr_elements = attr_to_elements[0];
+
+         Array<int> ess_tdofs;
+         Array<int> boundary_dofs;
+         fes.GetBoundaryLoopEdgeDofs(bdr_elements, ess_tdofs, boundary_dofs);
+
+         // Identify DOFs by global true DOF number, which is agreed upon by all
+         // ranks sharing the DOF, then compare the set selected anywhere with
+         // the set actually emitted in ess_tdof_list.
+         std::set<HYPRE_BigInt> selected, emitted;
+         for (int dof : boundary_dofs)
+         {
+            selected.insert(fes.GetGlobalTDofNumber(dof));
+         }
+         for (int i = 0; i < ess_tdofs.Size(); i++)
+         {
+            emitted.insert(fes.GetMyTDofOffset() + ess_tdofs[i]);
+         }
+
+         auto all_gather = [nranks](const std::set<HYPRE_BigInt> &s)
+         {
+            std::vector<HYPRE_BigInt> local(s.begin(), s.end());
+            int n = static_cast<int>(local.size()), total = 0;
+            std::vector<int> counts(nranks), bytes(nranks), displs(nranks);
+            MPI_Allgather(&n, 1, MPI_INT, counts.data(), 1, MPI_INT,
+                          MPI_COMM_WORLD);
+            constexpr int sz = sizeof(HYPRE_BigInt);
+            for (int r = 0; r < nranks; r++)
+            {
+               displs[r] = total * sz;
+               total += counts[r];
+               bytes[r] = counts[r] * sz;
+            }
+            std::vector<HYPRE_BigInt> all(total);
+            MPI_Allgatherv(local.data(), n * sz, MPI_BYTE, all.data(),
+                           bytes.data(), displs.data(), MPI_BYTE,
+                           MPI_COMM_WORLD);
+            return std::set<HYPRE_BigInt>(all.begin(), all.end());
+         };
+
+         // Gather both sets across all ranks. global_selected is every shared
+         // boundary DOF chosen on any rank; global_emitted is every true DOF
+         // actually placed in some rank's ess_tdof_list. A selected DOF missing
+         // from global_emitted is one that no rank owns and outputs, which is
+         // exactly the synchronization bug this test guards against.
+         const std::set<HYPRE_BigInt> global_selected = all_gather(selected);
+         const std::set<HYPRE_BigInt> global_emitted = all_gather(emitted);
+
+         int num_missing = 0;
+         for (auto gtdof : global_selected)
+         {
+            if (!global_emitted.count(gtdof)) { num_missing++; }
+         }
+
+         CAPTURE(orientation, nranks, global_selected.size(),
+                 global_emitted.size(), num_missing);
+         REQUIRE(num_missing == 0);
+      }
+   }
+}
+
+TEST_CASE("BoundaryEdgeDOFs2DLoopVertexDOFsPartitionInvariant",
+          "[Parallel][ParMesh][BoundaryEdgeDOFs]")
+{
+   // A closed boundary loop split between ranks must give the same result as
+   // the serial code. With a collection carrying vertex DOFs (ND_R2D), a vertex
+   // shared by two boundary segments is interior to the loop and must be
+   // dropped. When the two segments live on different ranks, each rank sees the
+   // vertex only once locally, so the occurrence parity has to be reconciled
+   // across the sharing group.
+   if (Mpi::WorldSize() < 2) { return; }
+
+   constexpr int order = 1;
+   ND_R2D_FECollection fec(order, 2);
+
+   // Serial reference result
+   Mesh serial_mesh = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false,
+                                            1.0, 1.0);
+   FiniteElementSpace serial_fes(&serial_mesh, &fec);
+
+   Array<int> serial_bdr_elements(serial_mesh.GetNBE());
+   for (int i = 0; i < serial_bdr_elements.Size(); i++)
+   {
+      serial_bdr_elements[i] = i;
+   }
+
+   Array<int> serial_boundary_dofs;
+
+   serial_fes.GetBoundaryLoopEdgeDofs(serial_bdr_elements, serial_boundary_dofs);
+
+   const int serial_count = serial_boundary_dofs.Size();
+
+   // Compare against several partitionings of the same mesh
+   const int num_procs = Mpi::WorldSize();
+   std::vector<std::vector<int>> partitionings;
+   {
+      Mesh probe = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false,
+                                         1.0, 1.0);
+      const int ne = probe.GetNE();
+
+      std::vector<int> block(ne), round_robin(ne);
+      for (int i = 0; i < ne; i++)
+      {
+         block[i] = (i < ne/2) ? 0 : num_procs-1;
+         round_robin[i] = i % num_procs;
+      }
+      partitionings.push_back(block);
+      partitionings.push_back(round_robin);
+   }
+
+   for (const auto &partition : partitionings)
+   {
+      Mesh mesh = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL, false,
+                                        1.0, 1.0);
+      ParMesh pmesh(MPI_COMM_WORLD, mesh, partition.data());
+      ParFiniteElementSpace pfes(&pmesh, &fec);
+
+      Array<int> local_bdr_elements(pmesh.GetNBE());
+      for (int i = 0; i < local_bdr_elements.Size(); i++)
+      {
+         local_bdr_elements[i] = i;
+      }
+
+      Array<int> ess_tdofs;
+      Array<int> local_boundary_dofs;
+      pfes.GetBoundaryLoopEdgeDofs(local_bdr_elements, ess_tdofs,
+                                   local_boundary_dofs);
+
+      // The true DOFs are owned by exactly one rank each, so summing the local
+      // counts gives a partition-independent global count.
+      int local_tdofs = ess_tdofs.Size();
+      int global_tdofs = 0;
+      MPI_Allreduce(&local_tdofs, &global_tdofs, 1, MPI_INT, MPI_SUM,
+                    MPI_COMM_WORLD);
+
+      CAPTURE(num_procs, serial_count, global_tdofs);
+      REQUIRE(global_tdofs == serial_count);
+   }
+}
+
+TEST_CASE("GroupCommunicatorReduceMarkedByGroupStride",
+          "[Parallel][GroupCommunicator]")
+{
+   // Regression test for the neighbor-major stride of the byGroup receive
+   // buffer: with more than one DOF in a group, the contributions to DOF i are
+   // at buf[j*nldofs + i], so reducing a single marked DOF must gather the
+   // strided values rather than reading a contiguous run.
+   const int rank = Mpi::WorldRank();
+   const int nranks = Mpi::WorldSize();
+   if (nranks < 3) { return; }
+
+   ListOfIntegerSets groups;
+
+   IntegerSet local_group(1);
+   local_group[0] = rank;
+   groups.Insert(local_group);
+
+   IntegerSet shared_group(nranks);
+   for (int r = 0; r < nranks; r++)
+   {
+      shared_group[r] = r;
+   }
+   groups.Insert(shared_group);
+
+   GroupTopology topology(MPI_COMM_WORLD);
+   topology.Create(groups, 4983);
+
+   GroupCommunicator comm(topology, GroupCommunicator::byGroup);
+   // Two DOFs in the same shared group, so the buffer stride is 2.
+   Array<int> ldof_group(2);
+   ldof_group = 1;
+   comm.Create(ldof_group);
+
+   Array<real_t> values(2);
+   values[0] = real_t(10.0) * rank + real_t(1.0);
+   values[1] = real_t(100.0) * rank + real_t(2.0);
+
+   Array<int> marker(2);
+   marker = 1;
+
+   comm.ReduceBegin(values.GetData());
+   comm.ReduceMarked<real_t>(values.GetData(), marker, 0,
+                             GroupCommunicator::Sum<real_t>);
+   comm.Bcast(values);
+
+   const real_t rank_sum = real_t(nranks) * real_t(nranks - 1) / real_t(2.0);
+   REQUIRE(values[0] == MFEM_Approx(real_t(10.0) * rank_sum + real_t(nranks)));
+   REQUIRE(values[1] == MFEM_Approx(real_t(100.0) * rank_sum +
+                                    real_t(2.0) * real_t(nranks)));
+}
+
+TEST_CASE("GroupCommunicatorMaxAbs", "[Parallel][GroupCommunicator]")
+{
+   const int rank = Mpi::WorldRank();
+   const int nranks = Mpi::WorldSize();
+   if (nranks < 2) { return; }
+
+   ListOfIntegerSets groups;
+
+   IntegerSet local_group(1);
+   local_group[0] = rank;
+   groups.Insert(local_group);
+
+   IntegerSet shared_group(nranks);
+   for (int r = 0; r < nranks; r++)
+   {
+      shared_group[r] = r;
+   }
+   groups.Insert(shared_group);
+
+   GroupTopology topology(MPI_COMM_WORLD);
+   topology.Create(groups, 4983);
+
+   GroupCommunicator comm(topology, GroupCommunicator::byGroup);
+   Array<int> ldof_group(2);
+   ldof_group = 1;
+   comm.Create(ldof_group);
+
+   // The group master (rank 0) reduces the peers' contributions into its own.
+   Array<real_t> values(2);
+   // DOF 0: equal magnitude across ranks with opposite signs, with the negative
+   // value held by the master, so the opposite-sign tie must still resolve
+   // deterministically to the positive value.
+   values[0] = (rank == 0) ? real_t(-5.0) : real_t(5.0);
+   // DOF 1: the largest magnitude is negative and held by a peer, so it must
+   // win over the master's smaller positive value and keep its sign.
+   values[1] = (rank == nranks - 1) ? real_t(-10.0) : real_t(5.0);
+
+   comm.Reduce<real_t>(values.GetData(), GroupCommunicator::MaxAbs<real_t>);
+   comm.Bcast(values);
+
+   REQUIRE(values[0] == MFEM_Approx(real_t(5.0)));
+   REQUIRE(values[1] == MFEM_Approx(real_t(-10.0)));
+}
+
+#endif // MFEM_USE_MPI
Author	SHA1	Message	Date
Tzanio Kolev	5581b0c079	Merge pull request #4983 from dnpham23/master Methods for extracting boundary edge DOFs from a 3D ParMesh and for synchronizing boundary conditions across processor boundaries	2026-08-20 15:23:36 -07:00
Dzung Pham	fcbbfd49b9	Address review feedback: Fix H1 DOF dedup and single-precision reductions	2026-08-19 17:16:02 -07:00
Tzanio Kolev	5b1cc58146	Merge branch 'master' into master	2026-08-19 16:57:41 -07:00
Dzung Pham	240b2d5461	Make GroupCommunicator::MaxAbs tie-breaking deterministic On an equal-magnitude tie, MaxAbs now keeps the more positive value, so opposite-sign ties resolve to the positive one regardless of accumulation order; a strictly larger magnitude still wins and keeps its sign. Update the doc comment, add a GroupCommunicatorMaxAbs unit test, and add a CHANGELOG entry.	2026-08-19 15:55:05 -07:00
Dzung Pham	962943a743	Address najlkin review: tidy API and includes - GetBoundaryElementsByAttribute: return std::vector<Array<int>> indexed to match bdr_attrs instead of an unordered_map keyed by attribute. - ParFiniteElementSpace::GetBoundaryLoopEdgeDofs: make ldof_marker an optional trailing output (Array<int>*, default nullptr). It is still computed internally and swapped into the caller's array only when requested; the same information is derivable from boundary_edge_dofs via ListToMarker(). - Remove unused <unordered_set>/<unordered_map> includes from the fespace/pfespace headers; include them directly in pfespace.cpp where they are actually used. - Update the boundary-edge-DOF unit tests for the new signatures.	2026-08-19 13:30:32 -07:00
Dzung Pham	67e329b0e0	Add missing newline at end of test_bdr_edgedof.cpp	2026-08-19 12:15:43 -07:00
Tzanio Kolev	51a60d6460	Merge branch 'master' into master	2026-08-19 09:27:22 -07:00
Dzung Pham	aa0a0f1cc7	Fix -Wshadow error in boundary edge DOF test helper The inner loop index in ComputeBoundaryLoopLength shadowed the outer loop variable introduced when the helper was switched to iterate over the Array<int> dof_edges; rename it to fix the -Werror,-Wshadow CI build.	2026-08-18 23:02:42 -07:00
Dzung Pham	60d4ee9bf5	Address review feedback: Array<int> interface, doc/style fixes - Redesign GetBoundaryLoopEdgeDofs/ComputeLoopEdgeOrientations to expose single-indexed Array<int> outputs instead of maps/sets (maps kept only as internal scratch); drop the vestigial dof_to_orientation output. - Move the integral-type check into the BitOR template via static_assert and remove the unused ReduceOp enum Reduce overload. - Remove redundant ParFiniteElementSpace overrides that duplicated the base implementations (mesh == pmesh for all Par constructors). - Document the boundary loop definition and MaxAbs tie-breaking semantics. - Add CHANGELOG entry; hardcode the expected loop-edge DOF count in the partition-invariance test.	2026-08-18 22:48:33 -07:00
Dzung Pham	fb2ee4d248	Merge remote-tracking branch 'upstream/master'	2026-08-05 22:53:12 -07:00
Dzung Pham	88b8ead3b9	Address review feedback Document that GroupCommunicator::MaxAbs returns the signed value with the largest absolute value, not the non-negative absolute value.	2026-08-05 18:19:17 -07:00
Dzung Pham	81b6d40fd1	Address review feedback - Mention float in the Sum/Min/Max/MaxAbs instantiation comments. - Wrap a long GetBdrElementFace call to stay within 80 columns. - Explain in GetBoundaryLoopEdgeDofs why occurrences of GetEdgeDofs are counted rather than collecting GetEdgeInteriorDofs: the count also resolves vertex DOFs (kept at loop endpoints, dropped when shared), which GetEdgeInteriorDofs omits entirely.	2026-08-05 16:06:39 -07:00
Dzung Pham	2256251e56	Address review feedback - Serial GetBoundaryLoopEdgeDofs (3D and 2D): keep an edge/vertex DOF only when it appears in exactly one selected boundary element, matching the parallel version and handling non-manifold junctions correctly. - Inline the GetBoundaryElementsByAttributeImpl and ComputeLoopEdgeOrientationsImpl free helpers into their member functions. - Make GroupCommunicator::ReduceOp an enum class and note that the reduction result is signed. - Rename DoF -> DOF throughout the PR's additions and reflow doc comments to 80 columns. - Remove the unused GeneratePartitionings helper; build test partitions in place; add explanatory comments to the partition-invariant and shared-DOF-ownership tests.	2026-08-05 14:58:00 -07:00
Dzung Pham	dae8a18c32	Address review feedback - communication.hpp: discard the BitOR_Op branch at compile time for non-integral types, which otherwise failed to link. - communication.cpp: fix the neighbor-major buffer stride in ReduceMarked() byGroup mode, which mixed values from different DoFs. - pfespace.cpp: reset ess_edge_list alongside ess_tdof_list so reuse does not leave stale entries. - pfespace.cpp: reconcile 2D boundary vertex DoF occurrence parity across sharing groups, so a loop split between ranks matches serial. - pfespace.cpp: synchronize the local marker before converting to true DoFs, so a selected shared DoF is not skipped by every rank. Adds three regression tests, each confirmed to fail without its fix.	2026-07-29 12:09:39 -07:00
Hugh Carson	753e02c1c8	Merge branch 'master' into master	2026-07-17 15:16:05 -04:00
Dzung Pham	8f5c7a0eca	Merge remote-tracking branch 'upstream/master'	2026-07-01 13:43:07 -07:00
Dzung Pham	b33de09ae7	clean up unused variables	2025-12-04 11:43:58 -05:00
Dzung Pham	cfdf4321e2	Apply code style formatting to boundary edge DoF files and tests	2025-12-04 01:25:33 -05:00
Dzung Pham	84055907cb	Merge branch 'master' into master	2025-12-03 22:56:34 -05:00
Dzung Pham	e7523c2f87	refactoring functions in fespace and pfespace	2025-09-15 22:36:21 -04:00
Dzung Pham	97bd40f2f0	extend GetBoundaryEdgeDofs to cover both 2D and 3D cases, rename it to GetBoundaryLoopEdgeDofs, remove SynchronizeBC, add unit test for 2D case	2025-09-13 13:25:17 -04:00
Dzung Pham	8f090cf9dc	update unit test to reduce # of cases, replace test mesh with a leaner one, add documentation and some minor changes	2025-09-11 01:20:49 -04:00
Dzung Pham	ad83ec9d69	Refactoring	2025-09-04 03:02:40 -04:00
Dzung Pham	1049100166	WIP - refactoring	2025-08-29 11:38:32 -07:00
Dzung Pham	157a2360b4	fix format, add missing documentation, fix variable shadow	2025-08-21 18:41:42 -07:00
Dzung Pham	adbbeaccb3	remove temporary test folder	2025-08-20 17:10:28 -04:00
Dzung Pham	65257818c5	add unit tests for edge dof extractions and synchronization	2025-08-20 17:10:28 -04:00
Dzung Pham	bedb6f4050	clean up and reorganize the tests for boundary edge dof methods	2025-08-20 17:10:28 -04:00
Dzung Pham	51ed5a1359	remove out-of-date manual patch files	2025-08-20 17:10:28 -04:00
Dzung Pham	2fcb7bd083	Update tests for boundary edge DoFs	2025-08-20 17:10:28 -04:00
Dzung Pham	7b336c9ea2	update cubetest	2025-08-20 17:10:28 -04:00
Dzung Pham	b7784c087e	generalize the edge dof extraction method to multiple boundary loops	2025-08-20 17:10:28 -04:00
Dzung Pham	92fe6ae2bc	WIP	2025-08-20 17:10:28 -04:00
Dzung Pham	8793fe8536	Add alternative method to synchronize only marked ldofs across processors.	2025-08-20 17:10:28 -04:00
Dzung Pham	9a20a3b1ec	Update unit tests and cubetest to use new GetBoundaryEdgeDoFs and ReduceMarked (if applicable)	2025-08-20 17:10:28 -04:00
Dzung Pham	472241ccbf	Modify GetBoundaryEdgeDoFs to return ldof markers, add new reduction method ReduceMarked	2025-08-20 17:10:28 -04:00
Dzung Pham	90c8967758	Clean up the cube test	2025-08-20 17:10:28 -04:00
Dzung Pham	4e7aaac06a	Add new synchronize method for applying loop BC on tdofs in parallel. Add new MaxAbs reduce operation	2025-08-20 17:10:28 -04:00
Dzung Pham	b3f9e66822	Update Makefile for new test, modify cube test to use edge dof extraction from mfem	2025-08-20 17:10:28 -04:00
Dzung Pham	ab36f48dee	Add method to compute perimeters of bounary loops	2025-08-20 17:10:28 -04:00
Dzung Pham	ad47e2ec05	Add new test for loop edge boundary condition	2025-08-20 17:10:28 -04:00
Dzung Pham	03c24ccddc	Add new mesh for testing	2025-08-20 17:10:28 -04:00
Dzung Pham	7754a02670	Add patch files to include edge dof methods into pfespace	2025-08-20 17:10:28 -04:00
Dzung Pham	7778f2840a	Add brute force unit tests	2025-08-20 17:10:28 -04:00
Dzung Pham	dbae2da272	Add methods to extract boundary edge dofs and loop orientations	2025-08-20 17:10:28 -04:00