diff --git a/_build_system/build_config.py b/_build_system/build_config.py
index 17595562..21e6d321 100644
--- a/_build_system/build_config.py
+++ b/_build_system/build_config.py
@@ -230,6 +230,9 @@ def process_cmd_options(command_obj, cfs):
             else:
                 setattr(command_obj, attr, False)
 
+    if len(cfs) != 0:
+       assert False, "Following build options are undefined: " +  ", ".join(list(cfs))
+
 
 def process_setup_options(command_obj, args):
     for item in args:
diff --git a/_build_system/build_pymfem.py b/_build_system/build_pymfem.py
index 140a37be..41406326 100644
--- a/_build_system/build_pymfem.py
+++ b/_build_system/build_pymfem.py
@@ -147,7 +147,19 @@ def check_new(ifile):
         return os.path.getmtime(ifile) > os.path.getmtime(wfile)
 
     def update_integrator_exts():
-        pwd = chdir(os.path.join(rootdir, 'mfem', 'common'))
+        commdir = os.path.join(rootdir, 'mfem', 'common')
+        pwd = chdir(commdir)
+        
+        for filename in ('bilininteg_ext.i', 'lininteg_ext.i'):
+            fid = open(filename, "w")
+            fid.close()
+            
+        os.chdir(pwd)
+        for filename in ['lininteg.i', 'bilininteg.i']:
+            command = [swig_command] + swigflag + serflag + [filename]
+            make_call(command)
+
+        os.chdir(commdir)
         command1 = [sys.executable, "generate_lininteg_ext.py"]
         command2 = [sys.executable, "generate_bilininteg_ext.py"]
         make_call(command1)
@@ -192,9 +204,10 @@ def update_header_exists(mfem_source):
         serflag.append('-I' + os.path.join(bglb.suitesparse_prefix,
                                            'include', 'suitesparse'))
 
-    for filename in ['lininteg.i', 'bilininteg.i']:
-        command = [swig_command] + swigflag + serflag + [filename]
-        make_call(command)
+    # generate xxinteg_ext.i
+    #    create empty files
+    #    run swig with *.i files
+    #    call update scripts
     update_integrator_exts()
 
     commands = []
@@ -208,6 +221,10 @@ def update_header_exists(mfem_source):
     with mp_pool:
         mp_pool.map(subprocess.run, commands)
 
+    #for c in commands:
+    #   print(" ".join(c))
+    #    subprocess.run(c)
+
     if not do_parallel:
         os.chdir(pwd)
         return
diff --git a/_build_system/build_utils.py b/_build_system/build_utils.py
index 0cff67ef..30821667 100644
--- a/_build_system/build_utils.py
+++ b/_build_system/build_utils.py
@@ -19,7 +19,7 @@
 __all__ = ["read_mfem_tplflags", "abspath", "external_install_prefix",
            "make_call", "chdir", "remove_files",
            "make", "make_install", "download", "gitclone",
-           "record_mfem_sha", "cmake", 
+           "record_mfem_sha", "cmake",
            "get_numpy_inc", "get_mpi4py_inc", "find_libpath_from_prefix",
            "clean_so", ]
 
@@ -102,18 +102,24 @@ def make_call(command, target='', force_verbose=False, env=None):
         env.update(os.environ)
 
     myverbose = bglb.verbose or force_verbose
-    if not myverbose:
-        kwargs['stdout'] = subprocess.DEVNULL
-        kwargs['stderr'] = subprocess.DEVNULL
+
+    kwargs['stdout'] = subprocess.PIPE
+    kwargs['stderr'] = subprocess.PIPE
 
     p = subprocess.Popen(command, **kwargs)
-    p.communicate()
+    stdout, stderr = p.communicate()
     if p.returncode != 0:
         if target == '':
             target = " ".join(command)
+
+        print("!!!!!!!!!!!!!!! make call faield !!!!!!!!!!!!!!!!!")
         print("Failed when calling command: " + target)
+        print(stdout)
+        print(stderr)
         raise subprocess.CalledProcessError(p.returncode,
                                             " ".join(command))
+    if myverbose:
+        print(stdout)
 
 
 def chdir(path):
@@ -246,7 +252,6 @@ def cmake(path, **kwargs):
     make_call(command)
 
 
-
 def get_numpy_inc():
 
     python = sys.executable
@@ -307,6 +312,7 @@ def find_libpath_from_prefix(lib, prefix0):
 
     return ''
 
+
 def clean_so(all=None):
     python = sys.executable
     command = [python, "setup.py", "clean"]
@@ -330,6 +336,3 @@ def clean_so(all=None):
     make_call(command)
 
     chdir(pwd)
-
-
-
diff --git a/examples/ex5-hdg.py b/examples/ex5-hdg.py
new file mode 100644
index 00000000..4e7e3690
--- /dev/null
+++ b/examples/ex5-hdg.py
@@ -0,0 +1,464 @@
+'''
+   MFEM example 5
+
+   See c++ version in the MFEM library for more detail
+'''
+import os
+import time
+import mfem.ser as mfem
+
+from mfem.ser import intArray
+from os.path import expanduser, join, dirname
+import numpy as np
+from numpy import sin, cos, exp
+
+use_umfpack = hasattr(mfem, "UMFPackSolver")
+
+
+def run(order=1,
+        meshfile='',
+        nx=0,
+        ny=0,
+        dg=False,
+        brt=False,
+        td=False,
+        hb=False,
+        rd=False,
+        sn=False,
+        visualization=False,
+        device='cpu',
+        numba=False,
+        pa=False):
+
+    def pFunc_exact(x):
+        xi = float(x[0])
+        yi = float(x[1])
+        zi = 0.0
+        if len(x) == 3:
+            zi = x[2]
+        from numpy import sin, cos, exp
+        return exp(xi)*sin(yi)*cos(zi)
+
+    class uFunc_ex(mfem.VectorPyCoefficient):
+        def EvalValue(self, x):
+            xi = float(x[0])
+            yi = float(x[1])
+            zi = 0.0
+            if len(x) == 3:
+                zi = x[2]
+            ret = [- exp(xi)*sin(yi)*cos(zi),
+                   - exp(xi)*cos(yi)*cos(zi)]
+            if len(x) == 3:
+                ret.append(exp(xi)*sin(yi)*sin(zi))
+            return ret
+
+    class pFunc_ex(mfem.PyCoefficient):
+        def EvalValue(self, x):
+            return pFunc_exact(x)
+
+    class fFunc(mfem.VectorPyCoefficient):
+        def EvalValue(self, x):
+            if len(x) == 3:
+                return [0., 0., 0.]
+            else:
+                return [0., 0.]
+
+    class gFunc(mfem.PyCoefficient):
+        def EvalValue(self, x):
+            if len(x) == 3:
+                return -pFunc_exact(x)
+            return 0.
+
+    class f_natural(mfem.PyCoefficient):
+        def EvalValue(self, x):
+            return -pFunc_exact(x)
+
+    device = mfem.Device(device)
+    device.Print()
+
+    if meshfile != "":
+        mesh = mfem.Mesh(meshfile, 1, 1)
+    else:
+        if ny <= 0:
+            ny = nx
+        print(nx, ny)
+        mesh = mfem.Mesh.MakeCartesian2D(nx, ny, mfem.Element.QUADRILATERAL)
+
+    dim = mesh.Dimension()
+
+    # 4. Refine the mesh to increase the resolution. In this example we do
+    #    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
+    #    largest number that gives a final mesh with no more than 10,000
+    #    elements.
+
+    ref_levels = int(np.floor(np.log(10000./mesh.GetNE())/np.log(2.)/dim))
+    for x in range(ref_levels):
+        mesh.UniformRefinement()
+
+    # 5. Define a finite element space on the mesh. Here we use the
+    #    Raviart-Thomas finite elements of the specified order.
+    if dg:
+        R_coll = mfem.L2_FECollection(order, dim, mfem.BasisType.GaussLobatto)
+    elif brt:
+        R_coll = mfem.BrokenRT_FECollection(order, dim)
+        R_coll_dg = mfem.L2_FECollection(order+1, dim)
+    else:
+        R_coll = mfem.RT_FECollection(order, dim)
+
+    W_coll = mfem.L2_FECollection(order, dim)
+
+    R_space = mfem.FiniteElementSpace(mesh, R_coll, dim if dg else 1)
+    if brt:
+        R_space_hg = mfem.FiniteElementSpace(mesh, R_coll_dg, dim)
+    W_space = mfem.FiniteElementSpace(mesh, W_coll)
+
+    darcy = mfem.DarcyForm(R_space, W_space)
+
+    # 6. Define the BlockStructure of the problem, i.e. define the array of
+    #    offsets for each variable. The last component of the Array is the sum
+    #    of the dimensions of each block.
+
+    block_offsets = darcy.GetOffsets()
+    offsets = block_offsets.ToList()
+    print("***********************************************************")
+    print("dim(R) = " + str(offsets[1] - offsets[0]))
+    print("dim(W) = " + str(offsets[2] - offsets[1]))
+    print("dim(R+W) = " + str(offsets[-1]))
+    print("***********************************************************")
+
+    # 7. Define the coefficients, analytical solution, and rhs of the PDE.
+
+    kcoeff = mfem.ConstantCoefficient(1.0)  # acoustic resistance
+    ikcoeff = mfem.RatioCoefficient(1., kcoeff)  # inverse acoustic resistance
+
+    fcoeff = fFunc(dim)
+    fnatcoeff = f_natural()
+    gcoeff = gFunc()
+    ucoeff = uFunc_ex(dim)
+    pcoeff = pFunc_ex()
+
+    # 8. Allocate memory (x, rhs) for the analytical solution and the right hand
+    #    side.  Define the GridFunction u,p for the finite element solution and
+    #    linear forms fform and gform for the right hand side.  The data
+    #    allocated by x and rhs are passed as a reference to the grid functions
+    #    (u,p) and the linear forms (fform, gform).
+
+    mt = mfem.MemoryType(device.GetDeviceMemoryType())
+    x = mfem.BlockVector(block_offsets, mt)
+    rhs = mfem.BlockVector(block_offsets, mt)
+
+    fform = mfem.LinearForm()
+    fform.Update(R_space, rhs.GetBlock(0), 0)
+
+    if dg:
+        fform.AddDomainIntegrator(mfem.VectorDomainLFIntegrator(fcoeff))
+        fform.AddBdrFaceIntegrator(
+            mfem.VectorBoundaryFluxLFIntegrator(fnatcoeff))
+    else:
+        fform.AddDomainIntegrator(mfem.VectorFEDomainLFIntegrator(fcoeff))
+        if brt:
+            fform.AddBdrFaceIntegrator(
+                mfem.VectorFEBoundaryFluxLFIntegrator(fnatcoeff))
+        else:
+            fform.AddBoundaryIntegrator(
+                mfem.VectorFEBoundaryFluxLFIntegrator(fnatcoeff))
+
+    fform.Assemble()
+    fform.SyncAliasMemory(rhs)
+
+    gform = mfem.LinearForm()
+    gform.Update(W_space, rhs.GetBlock(1), 0)
+    gform.AddDomainIntegrator(mfem.DomainLFIntegrator(gcoeff))
+    gform.Assemble()
+    gform.SyncAliasMemory(rhs)
+
+    # 9. Assemble the finite element matrices for the Darcy operator
+    #                            D = [ M  B^T ]
+    #                                [ B   0  ]
+    #     where:
+    #
+    #     M = \int_\Omega k u_h \cdot v_h d\Omega   u_h, v_h \in R_h
+    #     B   = -\int_\Omega \div u_h q_h d\Omega   u_h \in R_h, q_h \in W_h
+    mVarf = darcy.GetFluxMassForm()
+    bVarf = darcy.GetFluxDivForm()
+
+    if dg:
+        mtVarf = darcy.GetPotentialMassForm()
+        mVarf.AddDomainIntegrator(mfem.VectorMassIntegrator(kcoeff))
+        bVarf.AddDomainIntegrator(mfem.VectorDivergenceIntegrator())
+        bVarf.AddInteriorFaceIntegrator(mfem.TransposeIntegrator(
+            mfem.DGNormalTraceIntegrator(-1.)))
+        mtVarf.AddInteriorFaceIntegrator(
+            mfem.HDGDiffusionIntegrator(ikcoeff, td))
+    else:
+        mVarf.AddDomainIntegrator(mfem.VectorFEMassIntegrator(kcoeff))
+        bVarf.AddDomainIntegrator(mfem.VectorFEDivergenceIntegrator())
+        if brt:
+            bVarf.AddInteriorFaceIntegrator(mfem.TransposeIntegrator(
+                mfem.DGNormalTraceIntegrator(-1.)))
+
+    # set hybridization / assembly level
+    chrono = time.time()
+
+    ess_flux_tdofs_list = mfem.intArray()
+    if hb:
+        trace_coll = mfem.DG_Interface_FECollection(order, dim)
+        trace_space = mfem.FiniteElementSpace(mesh, trace_coll)
+        darcy.EnableHybridization(trace_space,
+                                  mfem.NormalTraceJumpIntegrator(),
+                                  ess_flux_tdofs_list)
+    elif (rd and (dg or brt)):
+        darcy.EnableFluxReduction()
+
+    if pa:
+        mVarf.SetAssemblyLevel(mfem.AssemblyLevel_PARTIAL)
+
+    darcy.Assemble()
+
+    pDarcyOp = mfem.OperatorHandle()
+    X = mfem.Vector()
+    B = mfem.Vector()
+    x.Assign(0.0)
+    darcy.FormLinearSystem(ess_flux_tdofs_list, x, rhs, pDarcyOp, X, B)
+
+    print("Assembly took " + str(time.time() - chrono))
+
+    maxIter = 1000
+    rtol = 1e-6
+    atol = 1e-10
+
+    if hb or (rd and (dg or brt)):
+        # 10. Construct the preconditioner
+        prec = mfem.GSSmoother(pDarcyOp.AsSparseMatrix())
+
+        # 11. Solve the linear system with GMRES.
+        #     Check the norm of the unpreconditioned residual.
+        chrono = time.time()
+        solver = mfem.GMRESSolver()
+        solver.SetAbsTol(atol)
+        solver.SetRelTol(rtol)
+        solver.SetMaxIter(maxIter)
+        solver.SetOperator(pDarcyOp)
+        solver.SetPreconditioner(prec)
+        solver.SetPrintLevel(1)
+
+        solver.Mult(B, X)
+        darcy.RecoverFEMSolution(X, rhs, x)
+
+        if solver.GetConverged():
+            print("GMRES converged in " + str(solver.GetNumIterations()) +
+                  " iterations with a residual norm of " + str(solver.GetFinalNorm()))
+        else:
+            print("GMRES did not converge in " + str(solver.GetNumIterations()) +
+                  " iterations with a residual norm of " + str(solver.GetFinalNorm()))
+
+        print("GMRES solver took " + str(time.time() - chrono) + "s.\n")
+
+    else:
+        # 10. Construct the operators for preconditioner
+        #
+        #                 P = [ diag(M)         0         ]
+        #                     [  0       B diag(M)^-1 B^T ]
+        #
+        #     Here we use Symmetric Gauss-Seidel to approximate the inverse of the
+        #     pressure Schur Complement
+        MinvBt = mfem.SparseMatrix()
+        Md = mfem.Vector()
+        darcyPrec = mfem.BlockDiagonalPreconditioner(block_offsets)
+
+        if pa:
+            mVarf.AssembleDiagonal(Md)
+            Md_host = Md.HostRead()
+            intMd = mfem.Vector(mVarf.Height())
+
+            intMd.GetDataArray()[:] = 1./Md_host.GetDataArray()
+
+            BMBt_diag = mfem.Vector(bVarf.Height())
+            bVarf.AssembleDiagonal_ADAt(invMd, BMBt_diag)
+
+            ess_tdof_list = intArray()
+
+            invM = mfem.OperatorJacobiSmoother(Md, ess_tdof_list)
+            invS = mfem.OperatorJacobiSmoother(BMBt_diag, ess_tdof_list)
+        else:
+            M = mfem.SparseMatrix(mVarf.SpMat())
+            M.GetDiag(Md)
+            Md.HostReadWrite()
+
+            B = mfem.SparseMatrix(bVarf.SpMat())
+            MinvBt = mfem.Transpose(B)
+
+            for i in range(Md.Size()):
+                MinvBt.ScaleRow(i, 1./Md[i])
+
+            S = mfem.Mult(B, MinvBt)
+            if dg:
+                Mtm = mfem.SparseMatrix(mtVarf.SpMat())
+                Snew = mfem.Add(Mtm, S)
+                S = Snew
+            invM = mfem.DSmoother(M)
+
+            if use_umfpack:
+                invS = mfem.UMFPackSolver(S)
+            else:
+                invS = mfem.GSSmoother(S)
+
+    invM.iterative_mode = False
+    invS.iterative_mode = False
+
+    darcyPrec.SetDiagonalBlock(0, invM)
+    darcyPrec.SetDiagonalBlock(1, invS)
+
+    chrono = time.time()
+
+    solver = mfem.MINRESSolver()
+    solver.SetAbsTol(atol)
+    solver.SetRelTol(rtol)
+    solver.SetMaxIter(maxIter)
+    solver.SetOperator(pDarcyOp)
+    solver.SetPreconditioner(darcyPrec)
+    solver.SetPrintLevel(1)
+
+    solver.Mult(rhs, x)
+    solve_time = time.time() - chrono
+
+    if solver.GetConverged():
+        print("MINRES converged in " + str(solver.GetNumIterations()) +
+              " iterations with a residual norm of " + "{:g}".format(solver.GetFinalNorm()))
+    else:
+        print("MINRES did not converge in " + str(solver.GetNumIterations()) +
+              " iterations. Residual norm is " + "{:g}".format(solver.GetFinalNorm()))
+
+    print("MINRES solver took " + str(solve_time) + "s.")
+
+    # 12. Create the grid functions u and p. Compute the L2 error norms.
+
+    u = mfem.GridFunction()
+    p = mfem.GridFunction()
+    if dg:
+        u_broken = mfem.GridFunction(R_space, x.GetBlock(0), 0)
+        coeff = mfem.VectorGridFunctionCoefficient(u_broken)
+        u.SetSpace(R_space_dg)
+        u.ProjectCoefficient(coeff)
+    else:
+        u.MakeRef(R_space, x.GetBlock(0), 0)
+    p.MakeRef(W_space, x.GetBlock(1), 0)
+
+    order_quad = max(2, 2*order+1)
+
+    irs = [mfem.IntRules.Get(i, order_quad)
+           for i in range(mfem.Geometry.NumGeom)]
+
+    norm_p = mfem.ComputeLpNorm(2, pcoeff, mesh, irs)
+    norm_u = mfem.ComputeLpNorm(2, ucoeff, mesh, irs)
+    err_u = u.ComputeL2Error(ucoeff, irs)
+    err_p = p.ComputeL2Error(pcoeff, irs)
+
+    print("|| u_h - u_ex || / || u_ex || = " + "{:g}".format(err_u / norm_u))
+    print("|| p_h - p_ex || / || p_ex || = " + "{:g}".format(err_p / norm_p))
+
+    mesh.Print("ex5.mesh", 8)
+    u.Save("sol_u.gf", 8)
+    p.Save("sol_p.gf", 8)
+
+    visit_dc = mfem.VisItDataCollection("Example5", mesh)
+    visit_dc.RegisterField("velocity", u)
+    visit_dc.RegisterField("pressure", p)
+    visit_dc.Save()
+
+    paraview_dc = mfem.ParaViewDataCollection("Example5", mesh)
+    paraview_dc.SetPrefixPath("ParaView")
+    paraview_dc.SetLevelsOfDetail(order)
+    paraview_dc.SetCycle(0)
+    paraview_dc.SetDataFormat(mfem.VTKFormat_BINARY)
+    paraview_dc.SetHighOrderOutput(True)
+    paraview_dc.SetTime(0.0)
+    paraview_dc.RegisterField("velocity", u)
+    paraview_dc.RegisterField("pressure", p)
+    paraview_dc.Save()
+
+    if visualization:
+        sout_u = mfem.socketstream("localhost", 19916)
+        sout_u.precision(8)
+        sout_u << "solution\n" << mesh << u << "window_title 'Velocity'\n"
+        sout_u << "keys Rljvvvvvmmc"
+        sout_u.flush()
+
+        sout_p = mfem.socketstream("localhost", 19916)
+        sout_p.precision(8)
+        sout_p << "solution\n" << mesh << p << "window_title 'Pressure'\n"
+        sout_p << "keys Rljmmc"
+        sout_p.flush()
+
+
+if __name__ == "__main__":
+    from mfem.common.arg_parser import ArgParser
+
+    parser = ArgParser(description='Ex5 (Darcy Problem)')
+    parser.add_argument('-m', '--mesh',
+                        default="",
+                        action='store', type=str,
+                        help='Mesh file to use.')
+    parser.add_argument("-nx", "--ncells-x",
+                        action='store', type=int, default=5,
+                        help="Number of cells in x.")
+    parser.add_argument("-ny", "--ncells-y",
+                        action='store', type=int, default=5,
+                        help="Number of cells in y.")
+    parser.add_argument('-o', '--order',
+                        action='store', default=1, type=int,
+                        help="Finite element order (polynomial degree).")
+    parser.add_argument("-dg", "--discontinuous",
+                        action='store_true',
+                        help="Enable DG elements for fluxes.")
+    parser.add_argument("-brt", "--broken-RT",
+                        action='store_true',
+                        help="Enable broken RT elements for fluxes.")
+    parser.add_argument("-td", "--stab_diff",
+                        action='store', type=float, default=0.5,
+                        help="Diffusion stabilization factor (1/2=default)")
+    parser.add_argument("-hb", "--hybridization",
+                        action='store_true',
+                        help="Enable hybridization.")
+    parser.add_argument("-rd", "--reduction",
+                        action='store_true',
+                        help="Enable reduction of DG flux.")
+    parser.add_argument("-sn", "--solution-norm",
+                        action='store_true',
+                        help="Solution norm (0=native, 1=flux, 2=potential).")
+    parser.add_argument("-pa", "--partial-assembly",
+                        action='store_true',
+                        help="Enable Partial Assembly.")
+    parser.add_argument("-d", "--device",
+                        default="cpu", type=str,
+                        help="Device configuration string, see Device::Configure().")
+    parser.add_argument('-vis', '--visualization',
+                        action='store_true',
+                        help='Enable GLVis visualization')
+
+    args = parser.parse_args()
+    parser.print_options(args)
+
+    order = args.order
+    if args.mesh != "":
+        meshfile = expanduser(
+            join(dirname(__file__), '..', 'data', args.mesh))
+    else:
+        meshfile = ""
+    visualization = args.visualization
+    device = args.device
+    pa = args.partial_assembly
+
+    run(order=order,
+        meshfile=meshfile,
+        nx=args.ncells_x,
+        ny=args.ncells_y,
+        dg=args.discontinuous,
+        brt=args.broken_RT,
+        td=args.stab_diff,
+        hb=args.hybridization,
+        rd=args.reduction,
+        sn=args.solution_norm,
+        visualization=visualization,
+        device=device,
+        pa=pa)
diff --git a/mfem/_ser/bilininteg_hdg.i b/mfem/_ser/bilininteg_hdg.i
new file mode 100644
index 00000000..708de0d2
--- /dev/null
+++ b/mfem/_ser/bilininteg_hdg.i
@@ -0,0 +1,29 @@
+//
+// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved.
+//
+%module(package="mfem._ser") bilininteg_hdg
+
+%{
+#include "mfem.hpp"
+#include "numpy/arrayobject.h"
+#include "../common/pyoperator.hpp"
+#include "../common/pyintrules.hpp"
+#include "../common/pybilininteg.hpp"
+#include "../common/pycoefficient.hpp"
+#include "../common/pynonlininteg.hpp"
+%}
+
+%include "../common/existing_mfem_headers.i"
+#ifdef FILE_EXISTS_FEM_DARCY_BILININTEG_HDG
+
+%init %{
+import_array();
+%}
+%include "exception.i"
+%include "../common/exception.i"
+
+%import "bilininteg.i"
+
+%include "fem/darcy/bilininteg_hdg.hpp"
+
+#endif
diff --git a/mfem/_ser/darcyform.i b/mfem/_ser/darcyform.i
new file mode 100644
index 00000000..b0e33c2c
--- /dev/null
+++ b/mfem/_ser/darcyform.i
@@ -0,0 +1,34 @@
+//
+// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved.
+//
+%module(package="mfem._ser") darcyform
+
+%{
+#include "mfem.hpp"
+#include "numpy/arrayobject.h"
+#include "../common/pyoperator.hpp"
+#include "../common/pyintrules.hpp"
+#include "../common/pycoefficient.hpp"
+#include "../common/pybilininteg.hpp"  
+%}
+
+%include "../common/existing_mfem_headers.i"
+#ifdef FILE_EXISTS_FEM_DARCY_DARCYFORM
+
+%init %{
+import_array();
+%}
+%include "exception.i"
+%include "../common/exception.i"
+
+%import "fespace.i"
+%import "blockvector.i"
+%import "bilinearform.i"
+%import "nonlinearform.i"
+%import "darcyreduction.i"
+%import "darcyhybridization.i"
+
+%include "fem/darcy/darcyform.hpp"
+
+#endif
+
diff --git a/mfem/_ser/darcyhybridization.i b/mfem/_ser/darcyhybridization.i
new file mode 100644
index 00000000..8f628c9f
--- /dev/null
+++ b/mfem/_ser/darcyhybridization.i
@@ -0,0 +1,29 @@
+//
+// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved.
+//
+%module(package="mfem._ser") darcyhybridization
+
+%{
+#include "mfem.hpp"
+#include "numpy/arrayobject.h"
+#include "../common/pyoperator.hpp"
+#include "../common/pyintrules.hpp"
+#include "../common/pycoefficient.hpp"
+%}
+
+%include "../common/existing_mfem_headers.i"
+#ifdef FILE_EXISTS_FEM_DARCY_DARCYHYBRIDIZATION
+
+%init %{
+import_array();
+%}
+%include "exception.i"
+%include "../common/exception.i"
+
+ //%ignore mfem::DarcyHybridization::m_nlfi_u;
+%ignore mfem::DarcyHybridization;
+%import "nonlininteg.i"
+
+%include "fem/darcy/darcyhybridization.hpp"
+
+#endif
diff --git a/mfem/_ser/darcyreduction.i b/mfem/_ser/darcyreduction.i
new file mode 100644
index 00000000..52d77761
--- /dev/null
+++ b/mfem/_ser/darcyreduction.i
@@ -0,0 +1,33 @@
+//
+// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved.
+//
+%module(package="mfem._ser") darcyreduction
+
+%{
+#include "mfem.hpp"
+#include "numpy/arrayobject.h"
+#include "../common/pyoperator.hpp"
+#include "../common/pyintrules.hpp"
+#include "../common/pycoefficient.hpp"
+#include "../common/pybilininteg.hpp"
+#include "../common/pynonlininteg.hpp"
+%}
+
+%include "../common/existing_mfem_headers.i"
+#ifdef FILE_EXISTS_FEM_DARCY_DARCYREDUCTION
+
+%init %{
+import_array();
+%}
+%include "exception.i"
+%include "../common/exception.i"
+
+%import "fespace.i"
+%import "nonlinearform.i"
+%import "sparsemat.i"
+%import "blockmatrix.i"
+%import "blockvector.i"
+
+%include "fem/darcy/darcyreduction.hpp"
+
+#endif
diff --git a/mfem/_ser/estimators_hdg.i b/mfem/_ser/estimators_hdg.i
new file mode 100644
index 00000000..f910245f
--- /dev/null
+++ b/mfem/_ser/estimators_hdg.i
@@ -0,0 +1,28 @@
+//
+// Copyright (c) 2020-2025, Princeton Plasma Physics Laboratory, All rights reserved.
+//
+%module(package="mfem._ser") estimators_hdg
+
+%{
+#include "mfem.hpp"
+#include "numpy/arrayobject.h"
+#include "../common/pyoperator.hpp"
+#include "../common/pyintrules.hpp"
+#include "../common/pybilininteg.hpp"
+#include "../common/pycoefficient.hpp"
+%}
+
+%include "../common/existing_mfem_headers.i"
+#ifdef FILE_EXISTS_FEM_DARCY_ESTIMATORS_HDG
+
+%init %{
+import_array();
+%}
+%include "exception.i"
+%include "../common/exception.i"
+
+%import "bilininteg.i"
+
+%include "fem/darcy/estimators_hdg.hpp"
+
+#endif
diff --git a/mfem/_ser/nonlininteg.i b/mfem/_ser/nonlininteg.i
index 0b746101..118322e3 100644
--- a/mfem/_ser/nonlininteg.i
+++ b/mfem/_ser/nonlininteg.i
@@ -3,6 +3,7 @@
 //
 %module(package="mfem._ser", directors="1") nonlininteg
 %{
+#include "numpy/arrayobject.h"    
 #include "mfem.hpp"
 #include "../common/pycoefficient.hpp"
 #include "../common/pyoperator.hpp"
@@ -10,11 +11,11 @@
 #include "../common/pynonlininteg.hpp"
 using namespace mfem;
 %}
-/*
+
 %init %{
 import_array();
 %}
-*/
+
 /*
 %exception {
     try { $action }
diff --git a/mfem/_ser/setup.py b/mfem/_ser/setup.py
index 86b0c81a..b3d0d23a 100644
--- a/mfem/_ser/setup.py
+++ b/mfem/_ser/setup.py
@@ -110,7 +110,9 @@ def get_extensions():
                "attribute_sets", "arrays_by_name",
                "hyperbolic",
                "complex_densemat", "complexstaticcond", "complexweakform",
-               "bounds", "integrator"]
+               "bounds", "integrator",
+               "darcyform","estimators_hdg", "bilininteg_hdg",
+               "darcyreduction", "darcyhybridization"]
 
     if add_cuda == '1':
         from setup_local import cudainc
diff --git a/mfem/common/bilininteg_ext.i b/mfem/common/bilininteg_ext.i
deleted file mode 100644
index 0929cdc1..00000000
--- a/mfem/common/bilininteg_ext.i
+++ /dev/null
@@ -1,317 +0,0 @@
-namespace mfem {
-%pythonappend BilinearFormIntegrator::BilinearFormIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend TransposeIntegrator::TransposeIntegrator %{
-    if own_bfi_ == 1:  bfi_.thisown = 0
-%}
-
-%pythonappend LumpedIntegrator::LumpedIntegrator %{
-    if own_bfi_ == 1:  bfi_.thisown = 0
-%}
-
-%pythonappend InverseIntegrator::InverseIntegrator %{
-    if own_integ == 1:  integ.thisown = 0
-%}
-
-%pythonappend SumIntegrator::SumIntegrator %{
-    self.own_integs = own_integs
-%}
-
-%pythonappend MixedScalarIntegrator::MixedScalarIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedVectorIntegrator::MixedVectorIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedScalarVectorIntegrator::MixedScalarVectorIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedScalarMassIntegrator::MixedScalarMassIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedVectorProductIntegrator::MixedVectorProductIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedScalarDerivativeIntegrator::MixedScalarDerivativeIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedScalarWeakDerivativeIntegrator::MixedScalarWeakDerivativeIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedScalarDivergenceIntegrator::MixedScalarDivergenceIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedVectorDivergenceIntegrator::MixedVectorDivergenceIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedScalarWeakGradientIntegrator::MixedScalarWeakGradientIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedScalarCurlIntegrator::MixedScalarCurlIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedScalarWeakCurlIntegrator::MixedScalarWeakCurlIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedVectorMassIntegrator::MixedVectorMassIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedCrossProductIntegrator::MixedCrossProductIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedDotProductIntegrator::MixedDotProductIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedWeakGradDotIntegrator::MixedWeakGradDotIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedWeakDivCrossIntegrator::MixedWeakDivCrossIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedGradGradIntegrator::MixedGradGradIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedCrossGradGradIntegrator::MixedCrossGradGradIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedCurlCurlIntegrator::MixedCurlCurlIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedCrossCurlCurlIntegrator::MixedCrossCurlCurlIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedCrossCurlGradIntegrator::MixedCrossCurlGradIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedCrossGradCurlIntegrator::MixedCrossGradCurlIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedWeakCurlCrossIntegrator::MixedWeakCurlCrossIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedScalarWeakCurlCrossIntegrator::MixedScalarWeakCurlCrossIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedCrossGradIntegrator::MixedCrossGradIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedCrossCurlIntegrator::MixedCrossCurlIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedScalarCrossCurlIntegrator::MixedScalarCrossCurlIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedScalarCrossGradIntegrator::MixedScalarCrossGradIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedScalarCrossProductIntegrator::MixedScalarCrossProductIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedScalarWeakCrossProductIntegrator::MixedScalarWeakCrossProductIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedDirectionalDerivativeIntegrator::MixedDirectionalDerivativeIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedGradDivIntegrator::MixedGradDivIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedDivGradIntegrator::MixedDivGradIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedScalarWeakDivergenceIntegrator::MixedScalarWeakDivergenceIntegrator %{
-    self._coeff = vq
-%}
-
-%pythonappend MixedVectorGradientIntegrator::MixedVectorGradientIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedVectorCurlIntegrator::MixedVectorCurlIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedVectorWeakCurlIntegrator::MixedVectorWeakCurlIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedVectorWeakDivergenceIntegrator::MixedVectorWeakDivergenceIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend GradientIntegrator::GradientIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend DiffusionIntegrator::DiffusionIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MassIntegrator::MassIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend BoundaryMassIntegrator::BoundaryMassIntegrator %{
-    self._coeff = q
-%}
-
-%pythonappend ConvectionIntegrator::ConvectionIntegrator %{
-    self._coeff = q
-%}
-
-%pythonappend ConservativeConvectionIntegrator::ConservativeConvectionIntegrator %{
-    self._coeff = q
-%}
-
-%pythonappend GroupConvectionIntegrator::GroupConvectionIntegrator %{
-    self._coeff = q
-%}
-
-%pythonappend VectorMassIntegrator::VectorMassIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend VectorFEDivergenceIntegrator::VectorFEDivergenceIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend VectorFEWeakDivergenceIntegrator::VectorFEWeakDivergenceIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend VectorFECurlIntegrator::VectorFECurlIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend VectorFEBoundaryFluxIntegrator::VectorFEBoundaryFluxIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend DerivativeIntegrator::DerivativeIntegrator %{
-    self._coeff = q
-%}
-
-%pythonappend CurlCurlIntegrator::CurlCurlIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend VectorCurlCurlIntegrator::VectorCurlCurlIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend MixedCurlIntegrator::MixedCurlIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend VectorFEMassIntegrator::VectorFEMassIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend VectorDivergenceIntegrator::VectorDivergenceIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend DivDivIntegrator::DivDivIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend VectorDiffusionIntegrator::VectorDiffusionIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend ElasticityIntegrator::ElasticityIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend ElasticityComponentIntegrator::ElasticityComponentIntegrator %{
-    self._coeff = parent_
-%}
-
-%pythonappend DGTraceIntegrator::DGTraceIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend NonconservativeDGTraceIntegrator::NonconservativeDGTraceIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend DGDiffusionIntegrator::DGDiffusionIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend DGDiffusionBR2Integrator::DGDiffusionBR2Integrator %{
-    self._coeff = args
-%}
-
-%pythonappend DGElasticityIntegrator::DGElasticityIntegrator %{
-    self._coeff = args
-%}
-
-%pythonappend ScalarProductInterpolator::ScalarProductInterpolator %{
-    self._coeff = sc
-%}
-
-%pythonappend ScalarVectorProductInterpolator::ScalarVectorProductInterpolator %{
-    self._coeff = sc
-%}
-
-%pythonappend VectorScalarProductInterpolator::VectorScalarProductInterpolator %{
-    self._coeff = vc
-%}
-
-%pythonappend ScalarCrossProductInterpolator::ScalarCrossProductInterpolator %{
-    self._coeff = vc
-%}
-
-%pythonappend VectorCrossProductInterpolator::VectorCrossProductInterpolator %{
-    self._coeff = vc
-%}
-
-%pythonappend VectorInnerProductInterpolator::VectorInnerProductInterpolator %{
-    self._coeff = vc
-%}
-
-%pythonappend PyBilinearFormIntegrator::PyBilinearFormIntegrator %{
-    self._ir=ir
-%}
-
-%pythonappend SumIntegrator::AddIntegrator %{
-   if self.own_integs == 1: integ.thisown = 0
-%}
-}
diff --git a/mfem/common/generate_bilininteg_ext.py b/mfem/common/generate_bilininteg_ext.py
index 0d02825d..7267f841 100644
--- a/mfem/common/generate_bilininteg_ext.py
+++ b/mfem/common/generate_bilininteg_ext.py
@@ -41,6 +41,8 @@
             pass
         elif line.find("(self, vdim_)") != -1:
             pass
+        elif line.find("(self, sign_=1.)") != -1:
+            pass
         else:
             print(cname)
             print(line)
diff --git a/mfem/common/lininteg_ext.i b/mfem/common/lininteg_ext.i
deleted file mode 100644
index 7e02196b..00000000
--- a/mfem/common/lininteg_ext.i
+++ /dev/null
@@ -1,78 +0,0 @@
-namespace mfem {
-%pythonappend LinearFormIntegrator::LinearFormIntegrator %{
-    self._coeff = args
-%}
-%pythonappend DeltaLFIntegrator::DeltaLFIntegrator %{
-    self._coeff = args
-%}
-%pythonappend DomainLFIntegrator::DomainLFIntegrator %{
-    self._coeff = args
-%}
-%pythonappend DomainLFGradIntegrator::DomainLFGradIntegrator %{
-    self._coeff = QF
-%}
-%pythonappend BoundaryLFIntegrator::BoundaryLFIntegrator %{
-    self._coeff = QG
-%}
-%pythonappend BoundaryNormalLFIntegrator::BoundaryNormalLFIntegrator %{
-    self._coeff = QG
-%}
-%pythonappend BoundaryTangentialLFIntegrator::BoundaryTangentialLFIntegrator %{
-    self._coeff = QG
-%}
-%pythonappend VectorDomainLFIntegrator::VectorDomainLFIntegrator %{
-    self._ir=ir
-    self._coeff = QF
-%}
-%pythonappend VectorDomainLFGradIntegrator::VectorDomainLFGradIntegrator %{
-    self._coeff = QF
-%}
-%pythonappend VectorBoundaryLFIntegrator::VectorBoundaryLFIntegrator %{
-    self._coeff = QG
-%}
-%pythonappend VectorFEDomainLFIntegrator::VectorFEDomainLFIntegrator %{
-    self._coeff = F
-%}
-%pythonappend VectorFEDomainLFCurlIntegrator::VectorFEDomainLFCurlIntegrator %{
-    self._coeff = F
-%}
-%pythonappend VectorFEDomainLFDivIntegrator::VectorFEDomainLFDivIntegrator %{
-    self._coeff = QF
-%}
-%pythonappend VectorBoundaryFluxLFIntegrator::VectorBoundaryFluxLFIntegrator %{
-    self._ir=ir
-    self._coeff = f
-%}
-%pythonappend VectorFEBoundaryFluxLFIntegrator::VectorFEBoundaryFluxLFIntegrator %{
-    self._coeff = args
-%}
-%pythonappend VectorFEBoundaryNormalLFIntegrator::VectorFEBoundaryNormalLFIntegrator %{
-    self._coeff = f
-%}
-%pythonappend VectorFEBoundaryTangentLFIntegrator::VectorFEBoundaryTangentLFIntegrator %{
-    self._coeff = QG
-%}
-%pythonappend BoundaryFlowIntegrator::BoundaryFlowIntegrator %{
-    self._coeff = args
-%}
-%pythonappend DGDirichletLFIntegrator::DGDirichletLFIntegrator %{
-    self._coeff = args
-%}
-%pythonappend DGElasticityDirichletLFIntegrator::DGElasticityDirichletLFIntegrator %{
-    self._coeff = uD_
-%}
-%pythonappend WhiteGaussianNoiseDomainLFIntegrator::WhiteGaussianNoiseDomainLFIntegrator %{
-    self._coeff = QG
-%}
-%pythonappend VectorQuadratureLFIntegrator::VectorQuadratureLFIntegrator %{
-    self._ir=ir
-    self._coeff = vqfc
-%}
-%pythonappend QuadratureLFIntegrator::QuadratureLFIntegrator %{
-    self._ir=ir
-    self._coeff = qfc
-%}
-%pythonappend PyLinearFormIntegrator::PyLinearFormIntegrator %{
-    self._ir=ir
-%}
-}
\ No newline at end of file
diff --git a/mfem/ser.py b/mfem/ser.py
index a66f05e0..5c8b6d9c 100644
--- a/mfem/ser.py
+++ b/mfem/ser.py
@@ -85,6 +85,7 @@
 from  mfem._ser.submesh import *
 from  mfem._ser.transfermap import *
 from  mfem._ser.hyperbolic import *
+from mfem._ser.darcyform import *
 
 import mfem._ser.array as array
 import mfem._ser.blockoperator as blockoperator
@@ -103,6 +104,7 @@
 import mfem._ser.sparsemat as sparsemat
 import mfem._ser.tmop_modules as tmop
 
+
 #
 # modules not a part of standard build
 #