Higher-Order Multiscale Finite Element Method in C++

A C++/Intel MKL implementation of higher-order Multiscale Finite Element Methods (MsFEM) for two-dimensional elliptic equations, including oversampling, oscillatory boundary conditions, and an optional Petrov–Galerkin formulation.

Problem

The code targets heterogeneous diffusion problems on the unit square:

$$-\nabla\cdot\left(\kappa_\varepsilon(x,y)\nabla u\right)=f \quad\text{in }(0,1)^2,$$

with Dirichlet boundary conditions. Six coefficient examples are implemented in src/get_kappa.cpp, including smooth and rapidly oscillatory media. The oscillation scale is currently ε = 1/32.

What “higher-order” and “oversampling” mean here

• Higher-order basis: each coarse element uses tensor-product polynomial nodes; the current driver uses order 3 and accepts orders 1–4.
• Multiscale basis: local fine-grid solves incorporate κε into the coarse basis instead of using coefficient-independent polynomial functions.
• Oversampling: each local basis problem is solved on a region enlarged beyond the target coarse cell, then restricted back to the cell. This reduces artificial boundary effects from the local solves.
• Oscillatory boundary option: local basis functions can use coefficient-aware oscillatory boundary data instead of only polynomial traces.
• Petrov–Galerkin option: the trial and test spaces can differ during coarse assembly.

Current default experiment

The parameters are defined near the top of main_Oversampling.cpp:

Parameter	Default	Meaning
finemesh	128	fine elements per axis
Nx, Ny	16, 16	coarse elements per axis
order	3	polynomial order
example	2	coefficient/right-hand-side case
osci	true	use oscillatory local boundary data
pg	false	use standard Galerkin assembly
over	ceil(nx/4)	oversampling width in fine cells
parallel_num	16	requested OpenMP threads for the cell loop

finemesh must be divisible by both coarse dimensions. The main driver currently uses compile-time experiment constants; edit this block to change a run.

Build

Prerequisites:

• CMake 3.18 or newer;
• a C++17 compiler;
• Intel oneAPI MKL with its CMake package available;
• OpenMP (optional but recommended).

After sourcing the oneAPI environment:

source /opt/intel/oneapi/setvars.sh
cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build -j

If MKL is installed in a custom prefix, pass its package location through CMAKE_PREFIX_PATH:

cmake -S . -B build \
  -DCMAKE_BUILD_TYPE=Release \
  -DCMAKE_PREFIX_PATH=/path/to/oneapi/mkl/latest

No user-specific MKL path is hard-coded in the build files.

Run

The default problem allocates large dense work arrays. Remove the shell stack limit before running on Linux:

ulimit -s unlimited
./build/higher_order_msfem

Set OMP_NUM_THREADS consistently with parallel_num and the CPU allocation, for example:

export OMP_NUM_THREADS=16
./build/higher_order_msfem

Output and validation

The program prints:

• time spent constructing oscillatory boundary data;
• total wall-clock time;
• relative L2 errors when an analytical or reference solution is enabled;
• comparisons between reconstructed MsFEM and standard-FEM solutions in supported cases.

Examples 1 and 4 define analytical solutions in src/real_u.cpp. Other cases use f = 1 or coefficient-specific forcing and require an external reference solution before an error comparison is meaningful.

Source layout

.
├── main_Oversampling.cpp     experiment driver and coarse assembly
├── include/head.h            public function declarations
└── src/
    ├── get_kappa.cpp         heterogeneous coefficient cases
    ├── get_f.cpp             matching/source terms
    ├── solve_subfem2dho.cpp  higher-order local FEM solves
    ├── get_over_domain.cpp   oversampling domains
    ├── truncate_R.cpp        restrict oversampled bases
    ├── fix_basis_function.cpp
    └── ...                   quadrature, assembly, norms, and utilities

The local workflow is:

coarse cell
  -> enlarged oversampling region
  -> fine-grid coefficient and forcing
  -> local higher-order FEM solve
  -> restrict basis to the target cell
  -> assemble the global coarse operator
  -> solve and reconstruct on the fine grid

Reproducible comparisons

When comparing standard FEM, MsFEM, higher-order bases, or different oversampling widths, record:

• ε, coefficient example, and right-hand side;
• fine/coarse mesh sizes and polynomial order;
• boundary basis (osci) and formulation (pg);
• oversampling width;
• thread count, MKL version, compiler, and build type;
• relative L2 error and total/basis-construction time.

Limitations

• Experiment parameters are still source constants rather than command-line options.
• The default case is memory intensive.
• Only selected examples provide an analytical solution.
• The repository does not yet declare a software license; public visibility alone does not grant reuse rights.