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
Merge pull request #5439 from mfem/weighted-lor-transfer-v2
2026-08-20 15:23:36 -07:00 · 2026-08-20 11:38:06 -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
23 changed files with 4555 additions and 278 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
@@ -105,6 +115,8 @@ GPU computing

 - Added partial assembly support for MixedScalarWeakCrossProductIntegrator.

+- Added partial assembly support for MixedVectorGradientIntegrator for H1->RT.
+
 - Added support for device partial assembly CurlInterpolator.
  This supports 2D and 3D variants:
    2D H1 (out-of-plane) to RT (in-plane)
@@ -140,6 +152,14 @@ Miscellaneous
  using the new method ApplyDofSigns() in class ParFiniteElementSpace: the
  method will return immediately if no sign flips are needed.

+- Added support for coefficient-weighted LOR transfer in
+  L2ProjectionGridTransfer. The transfer conserves the weighted mass, for
+  example when transferring velocity while conserving density-weighted momentum.
+  This is illustrated in the lor-transfer and plor-transfer miniapps.
+
+- Added support for saving DataCollection output on the node-local storage,
+  instead of requiring that the filesystem is shared among all the ranks.
+
 API changes
 -----------
 - Removed ProjectGrad from 2D RT elements. Users should use ProjectCurl instead.
@@ -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
@@ -2013,6 +2013,9 @@ private:
   const DofToQuad *mapsC;         ///< Not owned. DOF-to-quad map, closed.
   const GeometricFactors *geom;   ///< Not owned
   int dim, ne, dofs1D, quad1D;
+   int op_entries;
+   FiniteElement::DerivType test_fetype =
+      mfem::FiniteElement::NONE; ///< Derivative type of the vector test space.
 };

 /** Class for integrating the bilinear form $a(u,v) := (Q \mathrm{curl}(u), v)$ in 3D and
@@ -1055,7 +1055,8 @@ public:

 typedef VectorCoefficient DiagonalMatrixCoefficient;

-/// Base class for Matrix Coefficients that optionally depend on time and space.
+/** Base class for matrix-valued coefficients that optionally depend on time
+    and space. */
 class MatrixCoefficient
 {
 protected:
@@ -1102,6 +1103,9 @@ public:
   /// the quadrature points. The matrix will be transposed or not according to
   /// the boolean argument @a transpose.
   ///
+   /// The stored entries use the same row/column convention as `Eval()`,
+   /// unless `transpose == true`, in which case `K^T` is stored instead.
+   ///
   /// The @a vdim of the QuadratureFunction should be equal to the height times
   /// the width of the matrix.
   virtual void Project(QuadratureFunction &qf, bool transpose=false);
@@ -38,9 +38,24 @@ int DataCollection::create_directory(const std::string &dir_name,
   // create directories recursively
   const char path_delim = '/';
   std::string::size_type pos = 0;
-   int err_flag;
+   int err_flag = 0;
 #ifdef MFEM_USE_MPI
   const ParMesh *pmesh = dynamic_cast<const ParMesh*>(mesh);
+   // In addition to the global root, let the lowest rank on each shared-memory
+   // node create the directory too, so that node-local (non-shared) filesystems
+   // get it on every node rather than only where the global root lives. On a
+   // shared filesystem the extra mkdir() hits EEXIST and is tolerated below.
+   bool node_root = true;
+   if (pmesh)
+   {
+      MPI_Comm node_comm;
+      MPI_Comm_split_type(pmesh->GetComm(), MPI_COMM_TYPE_SHARED, myid,
+                          MPI_INFO_NULL, &node_comm);
+      int node_rank;
+      MPI_Comm_rank(node_comm, &node_rank);
+      node_root = (node_rank == 0);
+      MPI_Comm_free(&node_comm);
+   }
 #endif

   do
@@ -52,7 +67,7 @@ int DataCollection::create_directory(const std::string &dir_name,
      err_flag = mkdir(subdir.c_str(), 0777);
      err_flag = (err_flag && (errno != EEXIST)) ? 1 : 0;
 #else
-      if (myid == 0 || pmesh == NULL)
+      if (node_root || pmesh == NULL)
      {
         err_flag = mkdir(subdir.c_str(), 0777);
         err_flag = (err_flag && (errno != EEXIST)) ? 1 : 0;
@@ -64,7 +79,8 @@ int DataCollection::create_directory(const std::string &dir_name,
 #ifdef MFEM_USE_MPI
   if (pmesh)
   {
-      MPI_Bcast(&err_flag, 1, MPI_INT, 0, pmesh->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE, &err_flag, 1, MPI_INT, MPI_MAX,
+                    pmesh->GetComm());
   }
 #endif

@@ -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);

@@ -1710,6 +1710,9 @@ void MixedScalarWeakCrossProductIntegrator::AssemblePA(
   pa_data.SetSize(dim * nq * ne, Device::GetMemoryType());
   PAHcurlDotSetup2D(quad1D, ne, trial_map_integral, ir->GetWeights(),
                     geom->J, coeff, pa_data);
+   // Match the extra sign introduced by the legacy assembled path's
+   // MixedScalarWeakCrossProductIntegrator::CalcShape().
+   pa_data *= -1.0;
 }

 void MixedScalarWeakCrossProductIntegrator::AddMultPA(const Vector &x,
@@ -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;
@@ -231,9 +231,11 @@ const Operator &InterpolationGridTransfer::BackwardOperator()

 L2ProjectionGridTransfer::L2Projection::L2Projection(
   const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_,
+   CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
   MemoryType d_mt_)
   : Operator(fes_lor_.GetVSize(), fes_ho_.GetVSize()),
-     fes_ho(fes_ho_), fes_lor(fes_lor_), d_mt(d_mt_)
+     fes_ho(fes_ho_), fes_lor(fes_lor_), coeff_ho(coeff_ho_),
+     coeff_lor(coeff_lor_), d_mt(d_mt_)
 { }

 void L2ProjectionGridTransfer::L2Projection::BuildHo2Lor(
@@ -263,12 +265,13 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
   IntegrationPointTransformation& ip_tr,
   DenseMatrix& M_mixed_el) const
 {
-   int order = fe_lor.GetOrder() + fe_ho.GetOrder() + tr_lor->OrderW();
-   const IntegrationRule* ir = &IntRules.Get(geom, order);
+   int order = fe_lor.GetOrder() + fe_ho.GetOrder() + tr_lor->OrderW() +
+               coeff_ho.order;
+   const IntegrationRule &ir = IntRules.Get(geom, order);
   M_mixed_el = 0.0;
-   for (int i = 0; i < ir->GetNPoints(); i++)
+   for (int i = 0; i < ir.GetNPoints(); i++)
   {
-      const IntegrationPoint& ip_lor = ir->IntPoint(i);
+      const IntegrationPoint& ip_lor = ir.IntPoint(i);
      IntegrationPoint ip_ho;
      ip_tr.Transform(ip_lor, ip_ho);
      Vector shape_lor(fe_lor.GetDof());
@@ -284,23 +287,23 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
      {
         w *= tr_lor->Weight();
      }
+      if (coeff_ho)
+      {
+         w *= coeff_ho.coeff->Eval(*tr_ho, ip_ho);
+      }
      shape_lor *= w;
      AddMultVWt(shape_lor, shape_ho, M_mixed_el);
   }
 }

-void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
-   Geometry::Type geom, const FiniteElement& fe_ho,
-   const FiniteElement& fe_lor, ElementTransformation* el_tr,
-   IntegrationPointTransformation& ip_tr,
+void L2ProjectionGridTransfer::L2Projection::ElemMixedEvaluation(
+   Geometry::Type geom, const FiniteElement& fe_ho, const FiniteElement& fe_lor,
+   IntegrationPointTransformation& ip_tr, const IntegrationRule& ir,
   DenseMatrix& B_L, DenseMatrix& B_H) const
 {
-   int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
-   const IntegrationRule* ir = &IntRules.Get(geom, order);
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
+   for (int i = 0; i < ir.GetNPoints(); i++)
   {
-      const IntegrationPoint& ip_lor = ir->IntPoint(i);
+      const IntegrationPoint& ip_lor = ir.IntPoint(i);
      IntegrationPoint ip_ho;

      // maps integration point ip_lor -> ip_ho
@@ -320,7 +323,6 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
         B_H(i, j) = shape_ho(j);
      }
   }
-
 }

 void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
@@ -328,10 +330,11 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
   const FiniteElementSpace& fes_lor_ea,
   Vector &M_LH, MemoryType d_mt_)
 {
-   Mesh* mesh_ho = fes_ho_ea.GetMesh();
-   Mesh* mesh_lor = fes_lor_ea.GetMesh();
-   int nel_ho = mesh_ho->GetNE();
-   int nel_lor = mesh_lor->GetNE();
+   Mesh &mesh_ho = *fes_ho_ea.GetMesh();
+   Mesh &mesh_lor = *fes_lor_ea.GetMesh();
+
+   const int nel_ho = mesh_ho.GetNE();
+   const int nel_lor = mesh_lor.GetNE();

   if (nel_ho == 0)
   {
@@ -339,11 +342,11 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
      return;
   }

-   const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
+   const CoarseFineTransformations& cf_tr = mesh_lor.GetRefinementTransforms();

   int nref_max = 0;
   Array<Geometry::Type> geoms;
-   mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
+   mesh_ho.GetGeometries(mesh_ho.Dimension(), geoms);
   for (int ig = 0; ig < geoms.Size(); ++ig)
   {
      Geometry::Type geom = geoms[ig];
@@ -360,130 +363,226 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
   {
      // Assume all HO elements are LOR in the same way
      const int iho = 0;
-      {
-         Array<int> lor_els;
-         ho2lor.GetRow(iho, lor_els);
-         int nref = ho2lor.RowSize(iho);
-
-         Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
-         const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
-         const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
-
-         // Allocate space for DenseTensors
-         ElementTransformation *el_tr = fes_lor_ea.GetElementTransformation(0);
-         int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
-         const IntegrationRule* ir_ea = &IntRules.Get(geom, order);
-         int qPts = ir_ea->GetNPoints();
-
-         // Containers for the basis functions sampled at quadrature points
-         B_L.SetSize(qPts, fe_lor.GetDof(), nref, d_mt);
-         B_H.SetSize(qPts, fe_ho.GetDof(), nref, d_mt);
-         D.SetSize(qPts, nref, nel_ho, d_mt);
-
-         const GeometricFactors *geo_facts =
-            mesh_lor->GetGeometricFactors(*ir_ea, GeometricFactors::DETERMINANTS);
-
-         MFEM_ASSERT(nel_ho*nref == nel_lor, "we expect nel_ho*nref == nel_lor");
-
-         // Setup data at quadrature points
-         // TODO add support for user coefficient
-         const auto W = Reshape(ir_ea->GetWeights().Read(), qPts);
-         const auto J = Reshape(geo_facts->detJ.Read(), qPts, nel_lor);
-         const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
-
-         mfem::forall(qPts * nref * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
-         {
-            const int q    = tid % qPts;
-            const int iref = (tid / qPts) % nref;
-            const int iho  = (tid / (qPts * nref)) % nel_ho;
-
-            const int lo_el_id = iref + nref*iho;
-            const real_t detJ = J(q, lo_el_id);
-
-            d_D(q, iref, iho) = W(q) * detJ;
-
-         });
-
-         emb_tr.SetIdentityTransformation(geom);
-         const DenseTensor &pmats = cf_tr.point_matrices[geom];
-
-         // Collect the basis functions
-         for (int iref = 0; iref < nref; ++iref)
-         {
-            int ilor = lor_els[iref];
-            // Now assemble the block-row of the mixed mass matrix associated
-            // with integrating HO functions against LOR functions on the LOR
-            // sub-element.
-
-            // Create the transformation that embeds the fine low-order element
-            // within the coarse high-order element in reference space
-            emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
-
-            DenseMatrix &b_lo = B_L(ilor);
-            DenseMatrix &b_ho = B_H(ilor);
-
-            ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, b_lo, b_ho);
-
-         } // loop over subcells of ho element
-         // end of quadrature point setup
-      }
-
-   } // completed setup of basis function and quadrature point
-
-   // Assemble mixed mass matrix
-   {
-      int iho = 0;
      Array<int> lor_els;
      ho2lor.GetRow(iho, lor_els);
-      int nref = ho2lor.RowSize(iho);
+      const int nref = ho2lor.RowSize(iho);
+      MFEM_VERIFY(nel_ho*nref == nel_lor, "we expect nel_ho*nref == nel_lor");
+
+      Geometry::Type geom = mesh_ho.GetElementBaseGeometry(iho);
+
+      emb_tr.SetIdentityTransformation(geom);
+      const DenseTensor &pmats = cf_tr.point_matrices[geom];

      const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
      const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
-      const int ndof_ho = fe_ho.GetDof();
-      const int ndof_lor = fe_lor.GetDof();

-      const int qPts = D.SizeI();
+      // Allocate space for DenseTensors
+      ElementTransformation &el_tr = *mesh_lor.GetTypicalElementTransformation();
+      const int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr.OrderW()
+                        + coeff_ho.order;
+      const IntegrationRule &ir_ea = IntRules.Get(geom, order);
+      const int qPts = ir_ea.GetNPoints();

-      M_LH.SetSize(ndof_lor*ndof_ho*nref*nel_ho, d_mt);
+      // Containers for the basis functions sampled at quadrature points
+      B_L.SetSize(qPts, fe_lor.GetDof(), nref, d_mt);
+      B_H.SetSize(qPts, fe_ho.GetDof(), nref, d_mt);
+      D.SetSize(qPts, nref, nel_ho, d_mt);

-      // Rows x columns
-      // Recall MFEM is column major
-      // rows x columns is inverted  - matrix is ndof_lor x ndof_ho
-      auto v_M_LH = Reshape(M_LH.Write(), ndof_lor, ndof_ho, nref,
-                            nel_ho);
+      const GeometricFactors *geo_facts =
+         mesh_lor.GetGeometricFactors(ir_ea, GeometricFactors::DETERMINANTS);

-      const int fe_ho_ndof = fe_ho.GetDof();
-      const int fe_lor_ndof = fe_lor.GetDof();
+      Vector coeff_vec(qPts*nel_lor);
+      coeff_vec.UseDevice(true);

-      auto d_B_L = Reshape(B_L.Read(), qPts, fe_lor_ndof, nref);
-      auto d_B_H = Reshape(B_H.Read(), qPts, fe_ho_ndof, nref);
-      auto d_D   = Reshape(D.Read(), qPts, nref, nel_ho);
+      const int dim = mesh_ho.Dimension();
+      const int nq1d = (int)floor(pow(ir_ea.Size(), 1.0/dim) + 0.5);
+      const int nref_1d = (int)floor(pow(nref, 1.0/dim) + 0.5);

-      mfem::forall(fe_ho_ndof*nref*nel_ho, [=] MFEM_HOST_DEVICE (int idx)
+      if (!coeff_ho)
      {
-         const int bh = idx % fe_ho_ndof;
-         const int iref = (idx / fe_ho_ndof) % nref;
-         const int iho = idx / fe_ho_ndof / nref;
-         // (B_lo_dofs x Q) x (Q x B_ho_dofs)
-         for (int bl = 0; bl < fe_lor_ndof; ++bl)
+         coeff_vec = 1.0;
+      }
+      else if (UsesTensorBasis(fes_ho) &&
+               nq1d*nref_1d <= DeviceDofQuadLimits::Get().MAX_Q1D)
+      {
+         // Fast coefficient evaluation for tensor-product case. We create a
+         // "composite" quadrature rule in the high-order element that is the
+         // union of the quadrature rules within each of the low-order-refined
+         // subelements.
+         //
+         // NOTE: if the integration rule order is high and there are many LOR
+         // subelements, this can create a very big quadrature rule. That is
+         // why we need to check that we do not exceed MAX_Q1D. If we do, then
+         // we fall back on the slower "legacy" evaluation.
+
+         // Construct the composite rule as a tensor-product of the 1D LOR rule.
+         IntegrationRule ir_ho = [&]()
         {
-            real_t dot = 0.0;
-            for (int qi=0; qi<qPts; ++qi)
+            IntegrationRule ir_ho_1d(nq1d * nref_1d);
+            for (int iref = 0; iref < nref_1d; ++iref)
            {
-               dot += d_B_L(qi, bl, iref) *  d_D(qi, iref, iho) * d_B_H(qi, bh, iref);
+               const real_t a = pmats(cf_tr.embeddings[iref].matrix)(0,0);
+               const real_t b = pmats(cf_tr.embeddings[iref].matrix)(0,1);
+               for (int iq = 0; iq < nq1d; ++iq)
+               {
+                  ir_ho_1d[iq + iref*nq1d].x = a + ir_ea[iq].x*(b - a);
+               }
+            }
+            if (dim == 1) { return ir_ho_1d; }
+            else if (dim == 2) { return IntegrationRule(ir_ho_1d, ir_ho_1d); }
+            else { return IntegrationRule(ir_ho_1d, ir_ho_1d, ir_ho_1d); }
+         }();
+
+         // Project the high-order coefficient on the high-order composite rule.
+         QuadratureSpace qs(mesh_ho, ir_ho);
+         CoefficientVector coeff_vec_ho(*coeff_ho.coeff, qs);
+
+         // Permute the coefficient values to the expected LOR ordering.
+         const int nq_ho = ir_ho.Size();
+         const auto d_Q_ho = Reshape(coeff_vec_ho.Read(), nq_ho, nel_ho);
+         const auto d_Q = Reshape(coeff_vec.Write(), qPts, nel_lor);
+
+         mfem::forall(nq_ho * nel_ho, [=] MFEM_HOST_DEVICE (int ii)
+         {
+            const int e_ho = ii / nq_ho;
+            const int iq_ho = ii % nq_ho;
+
+            int iq_tensor = iq_ho;
+            int iq_lor = 0;
+            int iref = 0;
+            int iq_stride = 1;
+            int iref_stride = 1;
+            const int nq_ho_1d = nq1d*nref_1d;
+            for (int d = 0; d < dim; ++d)
+            {
+               const int iq_ho_1d = iq_tensor % nq_ho_1d;
+               iq_tensor /= nq_ho_1d;
+
+               iq_lor += (iq_ho_1d % nq1d)*iq_stride;
+               iref += (iq_ho_1d / nq1d)*iref_stride;
+               iq_stride *= nq1d;
+               iref_stride *= nref_1d;
+            }
+            const int e_lor = iref + e_ho*nref;
+
+            d_Q(iq_lor, e_lor) = d_Q_ho(iq_ho, e_ho);
+         });
+      }
+      else
+      {
+         // Legacy/fallback coefficient evaluation for non-tensor-product cases
+         // or when the number of quadrature points is too large for the device
+         // kernels.
+         IntegrationPoint ip_ho;
+         for (int e_ho = 0; e_ho < nel_ho; ++e_ho)
+         {
+            ElementTransformation &ho_tr = *mesh_ho.GetElementTransformation(e_ho);
+            for (int iref = 0; iref < nref; ++iref)
+            {
+               const int e_lor = iref + e_ho*nref;
+               emb_tr.SetPointMat(pmats(cf_tr.embeddings[e_lor].matrix));
+
+               for (int iq = 0; iq < qPts; ++iq)
+               {
+                  const IntegrationPoint &ip_lor = ir_ea[iq];
+                  ip_tr.Transform(ip_lor, ip_ho);
+                  ho_tr.SetIntPoint(&ip_ho);
+                  coeff_vec[iq + e_lor*qPts] = coeff_ho.coeff->Eval(ho_tr, ip_ho);
+               }
            }
-            // column major storage
-            v_M_LH(bl, bh, iref, iho) = dot;
         }
+      }
+
+      // Setup data at quadrature points
+      const auto W = Reshape(ir_ea.GetWeights().Read(), qPts);
+      const auto J = Reshape(geo_facts->detJ.Read(), qPts, nel_lor);
+      const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
+      const auto d_Q = Reshape(coeff_vec.Read(), qPts, nel_lor);
+
+      mfem::forall(qPts * nref * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
+      {
+         const int q    = tid % qPts;
+         const int iref = (tid / qPts) % nref;
+         const int iho  = (tid / (qPts * nref)) % nel_ho;
+
+         const int lo_el_id = iref + nref*iho;
+         const real_t detJ = J(q, lo_el_id);
+
+         d_D(q, iref, iho) = W(q) * d_Q(q, lo_el_id) * detJ;
      });
-   } // end of mixed assembly mass matrix
+
+      // Collect the basis functions
+      for (int iref = 0; iref < nref; ++iref)
+      {
+         int ilor = lor_els[iref];
+         // Now assemble the block-row of the mixed mass matrix associated
+         // with integrating HO functions against LOR functions on the LOR
+         // sub-element.
+
+         // Create the transformation that embeds the fine low-order element
+         // within the coarse high-order element in reference space
+         emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
+
+         DenseMatrix &b_lo = B_L(ilor);
+         DenseMatrix &b_ho = B_H(ilor);
+
+         ElemMixedEvaluation(geom, fe_ho, fe_lor, ip_tr, ir_ea, b_lo, b_ho);
+      } // loop over subcells of ho element
+      // end of quadrature point setup
+   } // completed setup of basis function and quadrature point
+
+   // Assemble mixed mass matrix
+   int iho = 0;
+   Array<int> lor_els;
+   ho2lor.GetRow(iho, lor_els);
+   int nref = ho2lor.RowSize(iho);
+
+   const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
+   const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
+   const int ndof_ho = fe_ho.GetDof();
+   const int ndof_lor = fe_lor.GetDof();
+
+   const int qPts = D.SizeI();
+
+   M_LH.SetSize(ndof_lor*ndof_ho*nref*nel_ho, d_mt);
+
+   // Rows x columns
+   // Recall MFEM is column major
+   // rows x columns is inverted  - matrix is ndof_lor x ndof_ho
+   auto v_M_LH = Reshape(M_LH.Write(), ndof_lor, ndof_ho, nref,
+                         nel_ho);
+
+   const int fe_ho_ndof = fe_ho.GetDof();
+   const int fe_lor_ndof = fe_lor.GetDof();
+
+   auto d_B_L = Reshape(B_L.Read(), qPts, fe_lor_ndof, nref);
+   auto d_B_H = Reshape(B_H.Read(), qPts, fe_ho_ndof, nref);
+   auto d_D   = Reshape(D.Read(), qPts, nref, nel_ho);
+
+   mfem::forall(fe_ho_ndof*nref*nel_ho, [=] MFEM_HOST_DEVICE (int idx)
+   {
+      const int bh = idx % fe_ho_ndof;
+      const int iref = (idx / fe_ho_ndof) % nref;
+      const int iho = idx / fe_ho_ndof / nref;
+      // (B_lo_dofs x Q) x (Q x B_ho_dofs)
+      for (int bl = 0; bl < fe_lor_ndof; ++bl)
+      {
+         real_t dot = 0.0;
+         for (int qi=0; qi<qPts; ++qi)
+         {
+            dot += d_B_L(qi, bl, iref) *  d_D(qi, iref, iho) * d_B_H(qi, bh, iref);
+         }
+         // column major storage
+         v_M_LH(bl, bh, iref, iho) = dot;
+      }
+   });
 }

 L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
 (const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_,
+ CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
 const bool use_ea_, MemoryType d_mt_)
-   : L2Projection(fes_ho_, fes_lor_, d_mt_),
-     use_ea(use_ea_)
+   : L2Projection(fes_ho_, fes_lor_, coeff_ho_, coeff_lor_, d_mt_), use_ea(use_ea_)
 {
   if (use_ea)
   {
@@ -559,7 +658,11 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
      DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
      DenseMatrix M_mixed(ndof_lor*nref, ndof_ho);

-      MassIntegrator mi;
+      MassIntegrator mi = [&]()
+      {
+         return coeff_lor ? MassIntegrator(*coeff_lor.coeff) : MassIntegrator();
+      }();
+
      DenseMatrix M_lor_el(ndof_lor, ndof_lor);
      DenseMatrixInverse Minv_lor_el(&M_lor_el);
      DenseMatrix M_lor(ndof_lor*nref, ndof_lor*nref);
@@ -577,6 +680,10 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
         // Assemble the low-order refined mass matrix and invert locally
         int ilor = lor_els[iref];
         ElementTransformation *tr_lor = fes_lor.GetElementTransformation(ilor);
+
+         const int order = 2*fe_lor.GetOrder() + tr_lor->OrderW() + coeff_lor.order;
+         mi.SetIntegrationRule(IntRules.Get(geom, order));
+
         mi.AssembleElementMatrix(fe_lor, *tr_lor, M_lor_el);
         M_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
         Minv_lor_el.Factor();
@@ -668,25 +775,22 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::EAL2ProjectionL2Space()
   // Need to compute M_L
   // Note: Using user-inputted M_LH IntegrationRule ir
   // (higher order than needed) in order to re-use coeff
-   MassIntegrator mi;
+   MassIntegrator mi = [&]()
+   {
+      return coeff_lor ? MassIntegrator(*coeff_lor.coeff) : MassIntegrator();
+   }();
+
+   const int order = 2*fes_lor.GetMaxElementOrder()
+                     + mesh_lor->GetTypicalElementTransformation()->OrderW()
+                     + coeff_lor.order;
+   mi.SetIntegrationRule(
+      IntRules.Get(mesh_lor->GetTypicalElementGeometry(), order));

   Vector M_ea_lor;
-   int ndof_lor;
-   int ndof_ho;
-   int nref;
-   {
-      int iho = 0;
-      Array<int> lor_els;
-      ho2lor.GetRow(iho, lor_els);
-      nref = ho2lor.RowSize(iho);
-
-      const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
-      const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
-      ndof_ho = fe_ho.GetDof();
-      ndof_lor = fe_lor.GetDof();
-
-      M_ea_lor.SetSize(ndof_lor*ndof_lor*nel_lor, d_mt);
-   }
+   const int ndof_lor = fes_lor.GetTypicalFE()->GetDof();
+   const int ndof_ho = fes_ho.GetTypicalFE()->GetDof();
+   const int nref = ho2lor.RowSize(0);
+   M_ea_lor.SetSize(ndof_lor*ndof_lor*nel_lor, d_mt);

   const bool add = false;
   mi.AssembleEA(fes_lor, M_ea_lor, add);
@@ -1032,8 +1136,9 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::EAProlongateTranspose(

 L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
   const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_,
+   CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
   const bool use_ea_, MemoryType d_mt_)
-   : L2Projection(fes_ho_, fes_lor_, d_mt_),
+   : L2Projection(fes_ho_, fes_lor_, coeff_ho_, coeff_lor_, d_mt_),
     use_ea(use_ea_)
 {

@@ -1092,8 +1197,9 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(

 L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
   const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor,
+   CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
   const bool use_ea_, MemoryType d_mt_)
-   : L2Projection(pfes_ho, pfes_lor, d_mt_),
+   : L2Projection(pfes_ho, pfes_lor, coeff_ho_, coeff_lor_, d_mt_),
     use_ea(use_ea_), pcg(pfes_ho.GetComm())
 {

@@ -1165,12 +1271,12 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetupPCG()

 void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
 {
-   Mesh* mesh_ho = fes_ho.GetMesh();
-   Mesh* mesh_lor = fes_lor.GetMesh();
-   int nel_ho = mesh_ho->GetNE();
-   int nel_lor = mesh_lor->GetNE();
-   int ndof_ho = fes_ho.GetNDofs();
-   int ndof_lor = fes_lor.GetNDofs();
+   Mesh &mesh_ho = *fes_ho.GetMesh();
+   Mesh &mesh_lor = *fes_lor.GetMesh();
+   const int nel_ho = mesh_ho.GetNE();
+   const int nel_lor = mesh_lor.GetNE();
+   const int ndof_ho = fes_ho.GetNDofs();
+   const int ndof_lor = fes_lor.GetNDofs();

   // If the local mesh is empty, skip all computations
   if (nel_ho == 0)
@@ -1178,11 +1284,11 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
      return;
   }

-   const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
+   const CoarseFineTransformations& cf_tr = mesh_lor.GetRefinementTransforms();

   int nref_max = 0;
   Array<Geometry::Type> geoms;
-   mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
+   mesh_ho.GetGeometries(mesh_ho.Dimension(), geoms);
   for (int ig = 0; ig < geoms.Size(); ++ig)
   {
      Geometry::Type geom = geoms[ig];
@@ -1205,7 +1311,8 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()

   BilinearForm Mho(fes_ho_scalar.get());
   Mho.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   Mho.AddDomainIntegrator(new MassIntegrator);
+   Mho.AddDomainIntegrator(coeff_ho ? new MassIntegrator(*coeff_ho.coeff)
+                           : new MassIntegrator);
   Mho.Assemble();

   // Processor local lumped Mass
@@ -1215,7 +1322,16 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()

   BilinearForm Mlor(fes_lor_scalar.get());
   Mlor.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   Mlor.AddDomainIntegrator(new MassIntegrator);
+   {
+      MassIntegrator *mi = coeff_lor ? new MassIntegrator(*coeff_lor.coeff)
+                           : new MassIntegrator;
+      const int order = 2*fes_lor.GetMaxElementOrder()
+                        + mesh_lor.GetTypicalElementTransformation()->OrderW()
+                        + coeff_lor.order;
+      mi->SetIntegrationRule(
+         IntRules.Get(mesh_lor.GetTypicalElementGeometry(), order));
+      Mlor.AddDomainIntegrator(mi);
+   }
   Mlor.Assemble();

   Vector ones_lor(Mlor.Width()); ones_lor = 1.0;
@@ -1228,15 +1344,14 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
   MixedMassEA(fes_ho, fes_lor, M_LH_ea, d_mt);

   // Set ownership
-   M_LH_local_op = new H1SpaceMixedMassOperator(fes_ho_scalar.get(),
-                                                fes_lor_scalar.get(),
-                                                &ho2lor,
-                                                &M_LH_ea);
+   M_LH.reset(new H1SpaceMixedMassOperator(fes_ho_scalar.get(),
+                                           fes_lor_scalar.get(),
+                                           &ho2lor,
+                                           &M_LH_ea));

   ML_inv_vea.reset(new H1SpaceLumpedMassOperator(fes_ho_scalar.get(),
                                                  fes_lor_scalar.get(),
                                                  ML_inv_ea));
-   M_LH.reset(M_LH_local_op);
   R.reset(new ProductOperator(ML_inv_vea.get(), M_LH.get(), false,
                               false));

@@ -1253,18 +1368,18 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
 void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
 (const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor)
 {
-   Mesh* mesh_ho = pfes_ho.GetParMesh();
-   Mesh* mesh_lor = pfes_lor.GetParMesh();
-   int nel_ho = mesh_ho->GetNE();
-   int nel_lor = mesh_lor->GetNE();
+   Mesh &mesh_ho = *pfes_ho.GetParMesh();
+   Mesh &mesh_lor = *pfes_lor.GetParMesh();
+   int nel_ho = mesh_ho.GetNE();
+   int nel_lor = mesh_lor.GetNE();
   int ndof_ho = pfes_ho.GetNDofs();
   int ndof_lor = pfes_lor.GetNDofs();

-   const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
+   const CoarseFineTransformations& cf_tr = mesh_lor.GetRefinementTransforms();

   int nref_max = 0;
   Array<Geometry::Type> geoms;
-   mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
+   mesh_ho.GetGeometries(mesh_ho.Dimension(), geoms);
   for (int ig = 0; ig < geoms.Size(); ++ig)
   {
      Geometry::Type geom = geoms[ig];
@@ -1287,7 +1402,8 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space

   ParBilinearForm pMho(pfes_ho_scalar.get());
   pMho.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   pMho.AddDomainIntegrator(new MassIntegrator);
+   pMho.AddDomainIntegrator(coeff_ho ? new MassIntegrator(*coeff_ho.coeff)
+                            : new MassIntegrator);
   pMho.Assemble();

   // Processor local lumped Mass
@@ -1297,7 +1413,16 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space

   ParBilinearForm pMlor(pfes_lor_scalar.get());
   pMlor.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   pMlor.AddDomainIntegrator(new MassIntegrator);
+   {
+      MassIntegrator *mi = coeff_lor ? new MassIntegrator(*coeff_lor.coeff)
+                           : new MassIntegrator;
+      const int order = 2*fes_lor.GetMaxElementOrder()
+                        + mesh_lor.GetTypicalElementTransformation()->OrderW()
+                        + coeff_lor.order;
+      mi->SetIntegrationRule(
+         IntRules.Get(mesh_lor.GetTypicalElementGeometry(), order));
+      pMlor.AddDomainIntegrator(mi);
+   }
   pMlor.Assemble();

   Vector ones_lor(pMlor.Width()); ones_lor = 1.0;
@@ -1570,7 +1695,7 @@ std::unique_ptr<SparseMatrix>>
         int ilor = lor_els[iref];
         ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);

-         int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
+         int order = 2 * fe_lor.GetOrder() + el_tr->OrderW() + coeff_lor.order;
         const IntegrationRule* ir = &IntRules.Get(geom, order);
         ML_el = 0.0;
         for (int i = 0; i < ir->GetNPoints(); ++i)
@@ -1578,7 +1703,13 @@ std::unique_ptr<SparseMatrix>>
            const IntegrationPoint& ip_lor = ir->IntPoint(i);
            fe_lor.CalcShape(ip_lor, shape_lor);
            el_tr->SetIntPoint(&ip_lor);
-            ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
+            real_t w = ip_lor.weight;
+            if (coeff_lor)
+            {
+               w *= coeff_lor.coeff->Eval(*el_tr, ip_lor);
+            }
+            shape_lor *= el_tr->Weight() * w;
+            ML_el += shape_lor;
         }
         fes_lor.GetElementDofs(ilor, dofs_lor);
         ML_inv.AddElementVector(dofs_lor, ML_el);
@@ -2024,8 +2155,8 @@ void L2ProjectionGridTransfer::BuildF()
   {
      if (!Parallel())
      {
-         F = new L2ProjectionH1Space(dom_fes, ran_fes,
-                                     use_ea, d_mt);
+         F = new L2ProjectionH1Space(
+            dom_fes, ran_fes, coeff_ho, coeff_lor, use_ea, d_mt);
      }
      else
      {
@@ -2034,15 +2165,15 @@ void L2ProjectionGridTransfer::BuildF()
            static_cast<mfem::ParFiniteElementSpace&>(dom_fes);
         const mfem::ParFiniteElementSpace& ran_pfes =
            static_cast<mfem::ParFiniteElementSpace&>(ran_fes);
-         F = new L2ProjectionH1Space(dom_pfes, ran_pfes,
-                                     use_ea, d_mt);
+         F = new L2ProjectionH1Space(
+            dom_pfes, ran_pfes, coeff_ho, coeff_lor, use_ea, d_mt);
 #endif
      }
   }
   else
   {
-      F = new L2ProjectionL2Space(dom_fes, ran_fes,
-                                  use_ea, d_mt);
+      F = new L2ProjectionL2Space(
+         dom_fes, ran_fes, coeff_ho, coeff_lor, use_ea, d_mt);
   }
 }

@@ -19,6 +19,8 @@
 #include "pfespace.hpp"
 #endif

+#include <cstddef>
+
 namespace mfem
 {

@@ -162,6 +164,18 @@ public:
 };


+struct CoefficientWithOrder
+{
+   Coefficient *coeff;
+   int order;
+   CoefficientWithOrder() : coeff(nullptr), order(0) { }
+   CoefficientWithOrder(std::nullptr_t) : coeff(nullptr), order(0) { }
+   CoefficientWithOrder(Coefficient &coeff_) : coeff(&coeff_), order(1) { }
+   CoefficientWithOrder(Coefficient &coeff_, int order_)
+      : coeff(&coeff_), order(order_) { }
+   operator bool() const { return coeff != nullptr; }
+};
+
 /** @brief Transfer data in L2 and H1 finite element spaces between a coarse
    mesh and an embedded refined mesh using L2 projection. */
 /** The forward, coarse-to-fine, transfer uses L2 projection. The backward,
@@ -207,6 +221,8 @@ public:
   protected:
      const FiniteElementSpace& fes_ho;
      const FiniteElementSpace& fes_lor;
+      CoefficientWithOrder coeff_ho;
+      CoefficientWithOrder coeff_lor;

      MemoryType d_mt;
      Array<int> offsets;
@@ -214,8 +230,15 @@ public:

      L2Projection(const FiniteElementSpace& fes_ho_,
                   const FiniteElementSpace& fes_lor_,
+                   CoefficientWithOrder coeff_ho_,
+                   CoefficientWithOrder coeff_lor_,
                   MemoryType d_mt_ = Device::GetHostMemoryType());

+      L2Projection(const FiniteElementSpace& fes_ho_,
+                   const FiniteElementSpace& fes_lor_,
+                   MemoryType d_mt_ = Device::GetHostMemoryType())
+         : L2Projection(fes_ho_, fes_lor_, nullptr, nullptr, d_mt_) { }
+
      void BuildHo2Lor(int nel_ho, int nel_lor,
                       const CoarseFineTransformations& cf_tr);

@@ -225,11 +248,11 @@ public:
                         IntegrationPointTransformation& ip_tr,
                         DenseMatrix& M_mixed_el) const;

-      void ElemMixedMass(Geometry::Type geom, const FiniteElement& fe_ho,
-                         const FiniteElement& fe_lor,
-                         ElementTransformation* el_tr,
-                         IntegrationPointTransformation& ip_tr,
-                         DenseMatrix& B_L, DenseMatrix& B_H) const;
+      void ElemMixedEvaluation(Geometry::Type geom, const FiniteElement& fe_ho,
+                               const FiniteElement& fe_lor,
+                               IntegrationPointTransformation& ip_tr,
+                               const IntegrationRule& ir,
+                               DenseMatrix& B_L, DenseMatrix& B_H) const;
   public:
      /* Returns the Mixed Mass M_LH via device element assembly by building the
      basis functions and data at the quadrature points. */
@@ -287,9 +310,17 @@ public:
   public:
      L2ProjectionL2Space(const FiniteElementSpace& fes_ho_,
                          const FiniteElementSpace& fes_lor_,
+                          CoefficientWithOrder coeff_ho_,
+                          CoefficientWithOrder coeff_lor_,
                          const bool use_ea_,
                          MemoryType d_mt_ = Device::GetHostMemoryType());

+      L2ProjectionL2Space(const FiniteElementSpace& fes_ho_,
+                          const FiniteElementSpace& fes_lor_,
+                          const bool use_ea_,
+                          MemoryType d_mt_ = Device::GetHostMemoryType())
+         : L2ProjectionL2Space(fes_ho_, fes_lor_, nullptr, nullptr, use_ea_, d_mt_) { }
+
      /*Same as above but assembles and stores R_ea, P_ea */
      void EAL2ProjectionL2Space();

@@ -356,13 +387,30 @@ public:
   public:
      L2ProjectionH1Space(const FiniteElementSpace &fes_ho_,
                          const FiniteElementSpace &fes_lor_,
+                          CoefficientWithOrder coeff_ho_,
+                          CoefficientWithOrder coeff_lor_,
                          const bool use_ea_,
                          MemoryType d_mt_ = Device::GetHostMemoryType());
+
+      L2ProjectionH1Space(const FiniteElementSpace& fes_ho_,
+                          const FiniteElementSpace& fes_lor_,
+                          const bool use_ea_,
+                          MemoryType d_mt_ = Device::GetHostMemoryType())
+         : L2ProjectionH1Space(fes_ho_, fes_lor_, nullptr, nullptr, use_ea_, d_mt_) { }
+
 #ifdef MFEM_USE_MPI
      L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_,
                          const ParFiniteElementSpace &pfes_lor_,
+                          CoefficientWithOrder coeff_ho_,
+                          CoefficientWithOrder coeff_lor_,
                          const bool use_ea_,
                          MemoryType d_mt_ = Device::GetHostMemoryType());
+
+      L2ProjectionH1Space(const ParFiniteElementSpace& fes_ho_,
+                          const ParFiniteElementSpace& fes_lor_,
+                          const bool use_ea_,
+                          MemoryType d_mt_ = Device::GetHostMemoryType())
+         : L2ProjectionH1Space(fes_ho_, fes_lor_, nullptr, nullptr, use_ea_, d_mt_) { }
 #endif
      /// Same as above but assembles action of R through 4 parts:
      ///   ( )  inv( lumped(M_L) ), which is a diagonal matrix (essentially a vector)
@@ -508,18 +556,38 @@ public:
      virtual ~L2Prolongation() { }
   };

+   /// Coefficient for the mixed L2 inner product.
+   CoefficientWithOrder coeff_ho;
+   /// Coefficient for the low-order L2 inner product.
+   CoefficientWithOrder coeff_lor;
   L2Projection   *F; ///< Forward, coarse-to-fine, operator
   L2Prolongation *B; ///< Backward, fine-to-coarse, operator
   bool force_l2_space;

 public:
+   /// Construct the unweighted L2 projection grid transfer.
   L2ProjectionGridTransfer(FiniteElementSpace &coarse_fes_,
                            FiniteElementSpace &fine_fes_,
                            bool force_l2_space_ = false,
                            MemoryType d_mt_ = Device::GetHostMemoryType()) // move to method
      : GridTransfer(coarse_fes_, fine_fes_),
-        F(NULL), B(NULL), force_l2_space(force_l2_space_)
-   { }
+        coeff_ho(nullptr), coeff_lor(nullptr), F(nullptr), B(nullptr),
+        force_l2_space(force_l2_space_) { }
+
+   /// @brief Construct the weighted L2 projection grid transfer.
+   ///
+   /// The low-order inner product is weighted by @a coeff_lor, and the mixed
+   /// inner product is weighted by @a coeff_ho.
+   L2ProjectionGridTransfer(FiniteElementSpace &coarse_fes_,
+                            FiniteElementSpace &fine_fes_,
+                            CoefficientWithOrder coeff_ho_,
+                            CoefficientWithOrder coeff_lor_,
+                            bool force_l2_space_ = false,
+                            MemoryType d_mt_ = Device::GetHostMemoryType()) // move to method
+      : GridTransfer(coarse_fes_, fine_fes_),
+        coeff_ho(coeff_ho_), coeff_lor(coeff_lor_), F(nullptr), B(nullptr),
+        force_l2_space(force_l2_space_) { }
+
   virtual ~L2ProjectionGridTransfer();

   const Operator &ForwardOperator() override;
@@ -527,6 +595,7 @@ public:
   const Operator &BackwardOperator() override;

   bool SupportsBackwardsOperator() const override;
+
 private:
   void BuildF();
 };
@@ -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;
@@ -33,6 +33,7 @@
 //
 // Sample runs:  lor-transfer
 //               lor-transfer -h1
+//               lor-transfer -ea -w
 //               lor-transfer -t
 //               lor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4
 //               lor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2
@@ -59,11 +60,12 @@ string direction;

 // Exact functions to project
 real_t RHO_exact(const Vector &x);
+real_t W_exact(const Vector &x);
+real_t weight(const Vector &x);

 // Helper functions
 void visualize(VisItDataCollection &, string, int, int, int visport = 19916);
-real_t compute_mass(FiniteElementSpace *, real_t, VisItDataCollection &,
-                    string);
+real_t compute_mass(GridFunction &, real_t, string, CoefficientWithOrder);

 int main(int argc, char *argv[])
 {
@@ -76,6 +78,7 @@ int main(int argc, char *argv[])
   bool useH1 = false;
   int visport = 19916;
   bool use_pointwise_transfer = false;
+   bool use_weighted_transfer = false;
   const char *device_config = "cpu";
   bool use_ea       = false;

@@ -98,6 +101,9 @@ int main(int argc, char *argv[])
   args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
                  "-no-t", "--dont-use-pointwise-transfer",
                  "Use pointwise transfer operators instead of L2 projection.");
+   args.AddOption(&use_weighted_transfer, "-w", "--use-weighted-transfer",
+                  "-no-w", "--dont-use-weighted-transfer",
+                  "Use coefficient-weighted L2 projection.");
   args.AddOption(&device_config, "-d", "--device",
                  "Device configuration string, see Device::Configure().");
   args.AddOption(&use_ea, "-ea", "--ea-version", "-no-ea",
@@ -107,6 +113,15 @@ int main(int argc, char *argv[])
   // Configure device
   Device device(device_config);

+   if (use_weighted_transfer && !use_pointwise_transfer)
+   {
+      if (problem != 5)
+      {
+         cout << "Switching to positive problem = 5 for weighted transfer.\n";
+      }
+      problem = 5;
+   }
+
   // Read the mesh from the given mesh file.
   Mesh mesh(mesh_file, 1, 1);
   int dim = mesh.Dimension();
@@ -138,6 +153,14 @@ int main(int argc, char *argv[])
   FiniteElementSpace fespace(&mesh, fec);
   FiniteElementSpace fespace_lor(&mesh_lor, fec_lor);

+   FunctionCoefficient weight_fn_coeff(weight);
+   CoefficientWithOrder weight_coeff;
+   if (use_weighted_transfer)
+   {
+      weight_coeff.coeff = &weight_fn_coeff;
+      weight_coeff.order = 2;
+   }
+
   GridFunction rho(&fespace);
   GridFunction rho_lor(&fespace_lor);

@@ -165,7 +188,7 @@ int main(int argc, char *argv[])
   rho.SetTrueVector();
   rho.SetFromTrueVector();

-   real_t ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO       ");
+   real_t ho_mass = compute_mass(rho, -1.0, "HO       ", weight_coeff);
   if (vis) { visualize(HO_dc, "HO", Wx, Wy, visport); Wx += offx; }

   GridTransfer *gt;
@@ -175,7 +198,8 @@ int main(int argc, char *argv[])
   }
   else
   {
-      gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
+      gt = new L2ProjectionGridTransfer(fespace, fespace_lor, weight_coeff,
+                                        weight_coeff);
   }

   // Configure element assembly for device acceleration
@@ -186,9 +210,44 @@ int main(int argc, char *argv[])
   // HO->LOR restriction
   direction = "HO -> LOR @ LOR";
   R.Mult(rho, rho_lor);
-   compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO)    ");
+   compute_mass(rho_lor, ho_mass, "R(HO)    ", weight_coeff);
   if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy, visport); Wx += offx; }

+   if (use_weighted_transfer && !use_pointwise_transfer)
+   {
+      // Transfer velocity while conserving rho-weighted momentum.
+      GridFunctionCoefficient rho_coeff(&rho);
+      GridFunctionCoefficient rho_lor_coeff(&rho_lor);
+      ProductCoefficient prod_coeff(weight_fn_coeff, rho_coeff);
+      ProductCoefficient prod_lor_coeff(weight_fn_coeff, rho_lor_coeff);
+      CoefficientWithOrder prod_weight(prod_coeff, order + 2);
+      CoefficientWithOrder prod_lor_weight(prod_lor_coeff, lorder + 2);
+
+      GridFunction w(&fespace), w_lor(&fespace_lor);
+      FunctionCoefficient W(W_exact);
+      w.ProjectCoefficient(W);
+
+      cout << '\n';
+      const real_t ho_momentum = compute_mass(w, -1.0, "rho w HO ", prod_weight);
+
+      L2ProjectionGridTransfer vel_gt(fespace, fespace_lor, prod_weight,
+                                      prod_lor_weight);
+      vel_gt.UseEA(use_ea);
+      vel_gt.ForwardOperator().Mult(w, w_lor);
+      compute_mass(w_lor, ho_momentum, "rho w LOR", prod_lor_weight);
+
+      if (vel_gt.SupportsBackwardsOperator())
+      {
+         GridFunction w_prev = w;
+         vel_gt.BackwardOperator().Mult(w_lor, w);
+         compute_mass(w, ho_momentum, "P(rho w) ", prod_weight);
+
+         w_prev -= w;
+         cout.precision(12);
+         cout << "|w - P(R(w))|_∞     = " << w_prev.Normlinf() << "\n\n";
+      }
+   }
+
   if (gt->SupportsBackwardsOperator())
   {
      const Operator &P = gt->BackwardOperator();
@@ -196,7 +255,7 @@ int main(int argc, char *argv[])
      direction = "HO -> LOR @ HO";
      GridFunction rho_prev = rho;
      P.Mult(rho_lor, rho);
-      compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
+      compute_mass(rho, ho_mass, "P(R(HO)) ", weight_coeff);
      if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy, visport); Wx = 0; Wy += offy; }

      rho_prev -= rho;
@@ -218,7 +277,7 @@ int main(int argc, char *argv[])
   direction = "LOR -> HO @ LOR";
   rho_lor.ProjectCoefficient(RHO);
   GridFunction rho_lor_prev = rho_lor;
-   real_t lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR      ");
+   real_t lor_mass = compute_mass(rho_lor, -1.0, "LOR      ", weight_coeff);
   if (vis) { visualize(LOR_dc, "LOR", Wx, Wy, visport); Wx += offx; }

   if (gt->SupportsBackwardsOperator())
@@ -227,14 +286,14 @@ int main(int argc, char *argv[])
      // Prolongate to HO space
      direction = "LOR -> HO @ HO";
      P.Mult(rho_lor, rho);
-      compute_mass(&fespace, lor_mass, HO_dc, "P(LOR)   ");
+      compute_mass(rho, lor_mass, "P(LOR)   ", weight_coeff);
      if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy, visport); Wx += offx; }

      // Restrict back to LOR space. This won't give the original function because
      // the rho_lor doesn't necessarily live in the range of R.
      direction = "LOR -> HO @ LOR";
      R.Mult(rho, rho_lor);
-      compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
+      compute_mass(rho_lor, lor_mass, "R(P(LOR))", weight_coeff);
      if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy, visport); }

      rho_lor_prev -= rho_lor;
@@ -270,12 +329,26 @@ real_t RHO_exact(const Vector &x)
         return M_PI/2-atan(5*(2*x.Norml2()-1));
      case 4: // basis function
         return (x.Norml2() < 0.1) ? 1 : 0;
+      case 5: // positive function
+         return 2.0 + 2*x(0)*x(0) + 3*x(1)*x(1) - x(0)*x(1) + 0.1*sin(x.Norml2());
      default:
         return 1.0;
   }
 }


+real_t W_exact(const Vector &x)
+{
+   return x(1) + 0.25*cos(2*M_PI*x.Norml2());
+}
+
+
+real_t weight(const Vector &x)
+{
+   return x(0)*x(0) + x(1)*x(1) + 1.0;
+}
+
+
 void visualize(VisItDataCollection &dc, string prefix, int x, int y,
               int visport)
 {
@@ -292,21 +365,32 @@ void visualize(VisItDataCollection &dc, string prefix, int x, int y,
 }


-real_t compute_mass(FiniteElementSpace *L2, real_t massL2,
-                    VisItDataCollection &dc, string prefix)
+real_t compute_mass(GridFunction &gf, real_t oldmass, string prefix,
+                    CoefficientWithOrder mass_coeff)
 {
+   FiniteElementSpace &fes = *gf.FESpace();
+   Mesh &mesh = *fes.GetMesh();
+
+   // Integration order is a * (element order) + b.
+   const int a = 2;
+   const int b = mesh.GetTypicalElementTransformation()->OrderW() +
+                 mass_coeff.order;
+
   ConstantCoefficient one(1.0);
-   LinearForm lf(L2);
-   lf.AddDomainIntegrator(new DomainLFIntegrator(one));
+   Coefficient &coeff = mass_coeff ? *mass_coeff.coeff : one;
+   DomainLFIntegrator *integ = new DomainLFIntegrator(coeff, a, b);
+
+   LinearForm lf(&fes);
+   lf.AddDomainIntegrator(integ);
   lf.Assemble();

-   real_t newmass = lf(*dc.GetField("density"));
+   const real_t newmass = lf(gf);
   cout.precision(18);
   cout << space << " " << prefix << " mass   = " << newmass;
-   if (massL2 >= 0)
+   if (oldmass >= 0)
   {
      cout.precision(4);
-      cout << " ("  << fabs(newmass-massL2)*100/massL2 << "%)";
+      cout << " ("  << fabs(newmass-oldmass)*100/oldmass << "%)";
   }
   cout << endl;
   return newmass;
@@ -33,6 +33,7 @@
 //
 // Sample runs:  plor-transfer
 //               plor-transfer -h1
+//               plor-transfer -ea -w
 //               plor-transfer -t
 //               plor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4
 //               plor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2
@@ -59,11 +60,12 @@ string direction;

 // Exact functions to project
 real_t RHO_exact(const Vector &x);
+real_t W_exact(const Vector &x);
+real_t weight(const Vector &x);

 // Helper functions
 void visualize(VisItDataCollection &, string, int, int, int /* visport */);
-real_t compute_mass(ParFiniteElementSpace *, real_t, VisItDataCollection &,
-                    string);
+real_t compute_mass(ParGridFunction &, real_t, string, CoefficientWithOrder);

 int main(int argc, char *argv[])
 {
@@ -80,6 +82,7 @@ int main(int argc, char *argv[])
   bool useH1 = false;
   int visport = 19916;
   bool use_pointwise_transfer = false;
+   bool use_weighted_transfer = false;
   const char *device_config = "cpu";
   bool use_ea       = false;

@@ -102,6 +105,9 @@ int main(int argc, char *argv[])
   args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
                  "-no-t", "--dont-use-pointwise-transfer",
                  "Use pointwise transfer operators instead of L2 projection.");
+   args.AddOption(&use_weighted_transfer, "-w", "--use-weighted-transfer",
+                  "-no-w", "--dont-use-weighted-transfer",
+                  "Use coefficient-weighted L2 projection.");
   args.AddOption(&device_config, "-d", "--device",
                  "Device configuration string, see Device::Configure().");
   args.AddOption(&use_ea, "-ea", "--ea-version", "-no-ea",
@@ -112,6 +118,15 @@ int main(int argc, char *argv[])
   Device device(device_config);
   if (Mpi::Root()) { device.Print(); }

+   if (use_weighted_transfer && !use_pointwise_transfer)
+   {
+      if (problem != 5 && Mpi::Root())
+      {
+         cout << "Switching to positive problem = 5 for weighted transfer.\n";
+      }
+      problem = 5;
+   }
+
   // Read the mesh from the given mesh file.
   Mesh serial_mesh(mesh_file, 1, 1);
   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
@@ -154,6 +169,14 @@ int main(int argc, char *argv[])
   ParFiniteElementSpace fespace(&mesh, fec);
   ParFiniteElementSpace fespace_lor(&mesh_lor, fec_lor);

+   FunctionCoefficient weight_fn_coeff(weight);
+   CoefficientWithOrder weight_coeff;
+   if (use_weighted_transfer)
+   {
+      weight_coeff.coeff = &weight_fn_coeff;
+      weight_coeff.order = 2;
+   }
+
   ParGridFunction rho(&fespace);
   ParGridFunction rho_lor(&fespace_lor);

@@ -183,7 +206,7 @@ int main(int argc, char *argv[])
   rho.SetTrueVector();
   rho.SetFromTrueVector();

-   real_t ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO       ");
+   real_t ho_mass = compute_mass(rho, -1.0, "HO       ", weight_coeff);
   if (vis) { visualize(HO_dc, "HO", Wx, Wy, visport); Wx += offx; }

   GridTransfer *gt;
@@ -193,7 +216,8 @@ int main(int argc, char *argv[])
   }
   else
   {
-      gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
+      gt = new L2ProjectionGridTransfer(fespace, fespace_lor, weight_coeff,
+                                        weight_coeff);
   }

   // Configure element assembly for device acceleration
@@ -204,7 +228,7 @@ int main(int argc, char *argv[])
   // HO->LOR restriction
   direction = "HO -> LOR @ LOR";
   R.Mult(rho, rho_lor);
-   compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO)    ");
+   compute_mass(rho_lor, ho_mass, "R(HO)    ", weight_coeff);
   if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy, visport); Wx += offx; }
   auto global_max = [](const Vector& v)
   {
@@ -214,6 +238,47 @@ int main(int argc, char *argv[])
      return max;
   };

+   if (use_weighted_transfer && !use_pointwise_transfer)
+   {
+      // Transfer velocity while conserving rho-weighted momentum.
+      GridFunctionCoefficient rho_coeff(&rho);
+      GridFunctionCoefficient rho_lor_coeff(&rho_lor);
+      ProductCoefficient prod_coeff(weight_fn_coeff, rho_coeff);
+      ProductCoefficient prod_lor_coeff(weight_fn_coeff, rho_lor_coeff);
+      CoefficientWithOrder prod_weight(prod_coeff, order + 2);
+      CoefficientWithOrder prod_lor_weight(prod_lor_coeff, lorder + 2);
+
+      ParGridFunction w(&fespace), w_lor(&fespace_lor);
+      FunctionCoefficient W(W_exact);
+      w.ProjectCoefficient(W);
+
+      if (Mpi::Root()) { cout << '\n'; }
+      const real_t ho_momentum = compute_mass(w, -1.0, "rho w HO ", prod_weight);
+
+      L2ProjectionGridTransfer vel_gt(fespace, fespace_lor, prod_weight,
+                                      prod_lor_weight);
+      vel_gt.UseEA(use_ea);
+      vel_gt.ForwardOperator().Mult(w, w_lor);
+      compute_mass(w_lor, ho_momentum, "rho w LOR", prod_lor_weight);
+
+      if (vel_gt.SupportsBackwardsOperator())
+      {
+         ParGridFunction w_prev = w;
+         vel_gt.BackwardOperator().Mult(w_lor, w);
+         compute_mass(w, ho_momentum, "P(rho w) ", prod_weight);
+
+         w_prev -= w;
+         Vector w_prev_true(fespace.GetTrueVSize());
+         w_prev.GetTrueDofs(w_prev_true);
+         const real_t l_inf = global_max(w_prev_true);
+         if (Mpi::Root())
+         {
+            cout.precision(12);
+            cout << "|w - P(R(w))|_∞     = " << l_inf << "\n\n";
+         }
+      }
+   }
+
   if (gt->SupportsBackwardsOperator())
   {
      const Operator &P = gt->BackwardOperator();
@@ -221,7 +286,7 @@ int main(int argc, char *argv[])
      direction = "HO -> LOR @ HO";
      ParGridFunction rho_prev = rho;
      P.Mult(rho_lor, rho);
-      compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
+      compute_mass(rho, ho_mass, "P(R(HO)) ", weight_coeff);
      if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy, visport); Wx = 0; Wy += offy; }

      rho_prev -= rho;
@@ -263,7 +328,7 @@ int main(int argc, char *argv[])
   direction = "LOR -> HO @ LOR";
   rho_lor.ProjectCoefficient(RHO);
   ParGridFunction rho_lor_prev = rho_lor;
-   real_t lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR      ");
+   real_t lor_mass = compute_mass(rho_lor, -1.0, "LOR      ", weight_coeff);
   if (vis) { visualize(LOR_dc, "LOR", Wx, Wy, visport); Wx += offx; }

   if (gt->SupportsBackwardsOperator())
@@ -272,14 +337,14 @@ int main(int argc, char *argv[])
      // Prolongate to HO space
      direction = "LOR -> HO @ HO";
      P.Mult(rho_lor, rho);
-      compute_mass(&fespace, lor_mass, HO_dc, "P(LOR)   ");
+      compute_mass(rho, lor_mass, "P(LOR)   ", weight_coeff);
      if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy, visport); Wx += offx; }

      // Restrict back to LOR space. This won't give the original function because
      // the rho_lor doesn't necessarily live in the range of R.
      direction = "LOR -> HO @ LOR";
      R.Mult(rho, rho_lor);
-      compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
+      compute_mass(rho_lor, lor_mass, "R(P(LOR))", weight_coeff);
      if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy, visport); }

      rho_lor_prev -= rho_lor;
@@ -334,12 +399,26 @@ real_t RHO_exact(const Vector &x)
         return M_PI/2-atan(5*(2*x.Norml2()-1));
      case 4: // basis function
         return (x.Norml2() < 0.1) ? 1 : 0;
+      case 5: // positive function
+         return 2.0 + 2*x(0)*x(0) + 3*x(1)*x(1) - x(0)*x(1) + 0.1*sin(x.Norml2());
      default:
         return 1.0;
   }
 }


+real_t W_exact(const Vector &x)
+{
+   return x(1) + 0.25*cos(2*M_PI*x.Norml2());
+}
+
+
+real_t weight(const Vector &x)
+{
+   return x(0)*x(0) + x(1)*x(1) + 1.0;
+}
+
+
 void visualize(VisItDataCollection &dc, string prefix, int x, int y,
               int visport)
 {
@@ -358,23 +437,34 @@ void visualize(VisItDataCollection &dc, string prefix, int x, int y,
 }


-real_t compute_mass(ParFiniteElementSpace *L2, real_t massL2,
-                    VisItDataCollection &dc, string prefix)
+real_t compute_mass(ParGridFunction &gf, real_t oldmass, string prefix,
+                    CoefficientWithOrder mass_coeff)
 {
+   ParFiniteElementSpace &fes = *gf.ParFESpace();
+   Mesh &mesh = *fes.GetMesh();
+
+   // Integration order is a * (element order) + b.
+   const int a = 2;
+   const int b = mesh.GetTypicalElementTransformation()->OrderW() +
+                 mass_coeff.order;
+
   ConstantCoefficient one(1.0);
-   ParLinearForm lf(L2);
-   lf.AddDomainIntegrator(new DomainLFIntegrator(one));
+   Coefficient &coeff = mass_coeff ? *mass_coeff.coeff : one;
+   DomainLFIntegrator *integ = new DomainLFIntegrator(coeff, a, b);
+
+   ParLinearForm lf(&fes);
+   lf.AddDomainIntegrator(integ);
   lf.Assemble();

-   real_t newmass = lf(*dc.GetParField("density"));
+   const real_t newmass = lf(gf);
   if (Mpi::Root())
   {
      cout.precision(18);
      cout << space << " " << prefix << " mass   = " << newmass;
-      if (massL2 >= 0)
+      if (oldmass >= 0)
      {
         cout.precision(4);
-         cout << " ("  << fabs(newmass-massL2)*100/massL2 << "%)";
+         cout << " ("  << fabs(newmass-oldmass)*100/oldmass << "%)";
      }
      cout << endl;
   }
@@ -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
@@ -750,6 +750,89 @@ TEST_CASE("Hcurl/Hdiv Mixed PA Coefficient",
   }
 }

+TEST_CASE("Hcurl/Hdiv MixedVectorGradientPA",
+          "[GPU][PartialAssembly][Coefficient]")
+{
+   constexpr real_t tol = 4e-12;
+   dimension = GENERATE(2, 3);
+   // no coeff, scalar coeff, diagonal matrix coeff, full matrix coeff
+   auto coeffType = GENERATE(0, 1, 2, 3);
+   auto order = GENERATE(1, 2, 3);
+   // RT, ND
+   auto vFEType = GENERATE(0, 1);
+   CAPTURE(dimension, coeffType, order, vFEType);
+
+   const int ne = 3;
+   Mesh mesh = MakeCartesianNonaligned(dimension, ne);
+
+   H1_FECollection scalar_fec(order, dimension);
+   FiniteElementSpace s_fespace(&mesh, &scalar_fec);
+
+   std::unique_ptr<FiniteElementCollection> vector_fec;
+
+   switch (vFEType)
+   {
+      case 0:
+         vector_fec.reset(new RT_FECollection(order - 1, dimension));
+         break;
+      case 1:
+         vector_fec.reset(new ND_FECollection(order, dimension));
+         break;
+   }
+   FiniteElementSpace v_fespace(&mesh, vector_fec.get());
+
+   MixedBilinearForm pa_form(&s_fespace, &v_fespace);
+   pa_form.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   MixedBilinearForm fa_form(&s_fespace, &v_fespace);
+
+   std::unique_ptr<Coefficient> coeff;
+   std::unique_ptr<DiagonalMatrixCoefficient> dq_coeff;
+   std::unique_ptr<MatrixCoefficient> mq_coeff;
+   switch (coeffType)
+   {
+      case 0:
+         pa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator);
+         fa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator);
+         break;
+      case 1:
+         coeff.reset(new FunctionCoefficient(&coeffFunction));
+         pa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*coeff));
+         fa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*coeff));
+         break;
+      case 2:
+         dq_coeff.reset(new VectorFunctionCoefficient(dimension, &vectorCoeffFunction));
+         pa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*dq_coeff));
+         fa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*dq_coeff));
+         break;
+      case 3:
+         mq_coeff.reset(new MatrixFunctionCoefficient(
+                           dimension, &asymmetricMatrixCoeffFunction));
+         pa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*mq_coeff));
+         fa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*mq_coeff));
+         break;
+   }
+   pa_form.Assemble();
+   fa_form.Assemble();
+
+   GridFunction x_fa(&s_fespace), y_fa(&v_fespace), y_pa(&v_fespace);
+   x_fa.Randomize(1234);
+   REQUIRE(x_fa.Size() == pa_form.Width());
+   REQUIRE(x_fa.Size() == fa_form.Width());
+   REQUIRE(y_fa.Size() == fa_form.Height());
+   REQUIRE(y_pa.Size() == pa_form.Height());
+   pa_form.Mult(x_fa, y_pa);
+   fa_form.Mult(x_fa, y_fa);
+   y_pa -= y_fa;
+   REQUIRE(y_pa.Normlinf() <= tol);
+
+   GridFunction x_pa(&s_fespace);
+   y_fa.Randomize(1234);
+   pa_form.MultTranspose(y_fa, x_pa);
+   fa_form.MultTranspose(y_fa, x_fa);
+   x_pa -= x_fa;
+   REQUIRE(x_pa.Normlinf() <= tol);
+}
+
 TEST_CASE("3D Bilinear VectorFE Integrators PartialAssembly",
          "[BilinearFormIntegrator]"
          "[PartialAssembly]"
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
Tzanio Kolev	e0ef9a423c	Merge pull request #5439 from mfem/weighted-lor-transfer-v2 Add weighted LOR transfer	2026-08-20 11:38:06 -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
Tzanio Kolev	463cb07baf	Merge pull request #5456 from mfem/extra_gpu_em MixedVector gradient partial assembly	2026-08-19 13:19:13 -07:00
Andrew Ho	7e05f29325	changelog	2026-08-19 13:13:33 -07:00
Dzung Pham	67e329b0e0	Add missing newline at end of test_bdr_edgedof.cpp	2026-08-19 12:15:43 -07:00
Will Pazner	721d80b314	Fix LOR transfer miniapp integration on mixed meshes	2026-08-19 10:29:04 -07:00
Will Pazner	1bc33816f0	Merge remote-tracking branch 'origin/master' into weighted-lor-transfer-v2 # Conflicts: # CHANGELOG	2026-08-19 10:11:37 -07:00
Tzanio Kolev	51a60d6460	Merge branch 'master' into master	2026-08-19 09:27:22 -07:00
Tzanio Kolev	c661137756	Merge pull request #5415 from Sbozzolo/node-local-output-dirs Create node-local DataCollection output folders	2026-08-19 09:24:35 -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
Tzanio Kolev	5b3b486379	small fix	2026-08-18 18:30:11 -07:00
Andrew Ho	8789221a6b	review comments	2026-08-18 16:07:22 -07:00
Andrew Ho	69a7a605c0	Merge remote-tracking branch 'origin/extra_gpu_em' into extra_gpu_em	2026-08-18 15:42:24 -07:00
Andrew Ho	610ce458f6	Revert documentation comments	2026-08-18 15:41:41 -07:00
Tzanio Kolev	f7b6e0c0f0	Merge branch 'master' into weighted-lor-transfer-v2	2026-08-18 11:40:23 -07:00
Andrew Ho	7dea939ff8	Merge branch 'master' into extra_gpu_em	2026-08-18 10:57:46 -07:00
Andrew Ho	1ccb7bc613	Merge branch 'master' into extra_gpu_em	2026-08-18 10:57:24 -07:00
Andrew Ho	fbd217e8a4	Fixed bug for H1->RT	2026-08-18 10:35:32 -07:00
Andrew Ho	c11172b842	Added MultTranspose test It appears the bug for 3D H1->RT is tied to having NE > 1	2026-08-18 10:23:16 -07:00
Andrew Ho	0efbbd938a	Merge remote-tracking branch 'origin/extra_gpu_em' into extra_gpu_em	2026-08-18 08:27:01 -07:00
Andrew Ho	e5fae218af	Added PA tests for all coefficient types for MixedVectorGradientIntegrator Test seems to be failing for 3D RT	2026-08-18 00:32:39 -07:00
Andrew Ho	968dc0bfce	extra documentation from kris	2026-08-17 14:16:20 -07:00
Andrew Ho	8a88975532	Merge remote-tracking branch 'origin/master' into extra_gpu_em	2026-08-17 14:11:44 -07:00
Andrew Ho	b279e7f318	style	2026-08-17 14:11:30 -07:00
Andrew Ho	4b9d8b9247	Additional changes from Kris Beckwith	2026-08-17 13:48:19 -07:00
Tzanio Kolev	6f3ed5508a	Merge branch 'master' into node-local-output-dirs	2026-08-14 18:06:25 -07:00
Will Pazner	36be39433f	Add L2 projection transfer ctors without coefficients	2026-08-11 14:57:42 -07:00
Tzanio Kolev	56edc22b3a	Merge branch 'master' into weighted-lor-transfer-v2	2026-08-11 12:13:10 -07:00
Will Pazner	35b32b6a02	Make sure backwards operator is supported in plor-transfer	2026-08-11 11:36:18 -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
Will Pazner	ebbdd4bbb4	Merge remote-tracking branch 'origin/master' into weighted-lor-transfer-v2	2026-08-04 09:46:09 -07:00
Will Pazner	63c2be4ed6	Update CHANGELOG	2026-08-04 09:45:58 -07:00
Will Pazner	f324dd58d0	Add weighted LOR transfer sample runs	2026-08-04 09:42:21 -07:00
Will Pazner	626e4cc9c9	Check if backwards operator is supported in weighted LOR transfer	2026-08-04 09:42:12 -07:00
Will Pazner	f19dfabb75	Fix member variable shadowing	2026-08-04 09:37:10 -07:00
Will Pazner	b1b49cd3e9	Delete old comment	2026-08-04 09:35:33 -07:00
Will Pazner	2d33afe729	Remove unneeded ElementTransformation from ElemMixedEvaluation	2026-08-04 09:34:01 -07:00
Will Pazner	98b6f7c1cf	Remove unneeded gitignore	2026-08-04 09:31:36 -07:00
Will Pazner	7483034f7c	Add momentum-conserving weighted LOR transfer to miniapp	2026-08-04 09:30:04 -07:00
Will Pazner	48d16f7993	Remove standalone weighted LOR transfer miniapp	2026-08-03 16:09:27 -07:00
Will Pazner	1bfdf5bf31	Add option for weighted transfer in {lor,plor}_transfer miniapp	2026-08-03 16:09:06 -07:00
Will Pazner	5a28c20815	Add default constructor to CoefficientWithOrder	2026-08-03 16:09:06 -07:00
Will Pazner	3f38fc53f1	Add weighted LOR transfer example	2026-08-03 16:09:04 -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
Gabriele Bozzola	a2a14e8ad8	Add changelog entry	2026-07-28 12:12:23 -07:00
Gabriele Bozzola	98bbd8ad94	Merge branch 'master' into node-local-output-dirs	2026-07-28 12:10:45 -07:00
Will Pazner	2903d0f666	Support coefficient-weighted LOR transfer	2026-07-24 16:22:33 -07:00
Will Pazner	006e82f199	Remove unnecessary scope	2026-07-24 16:22:33 -07:00
Gabriele Bozzola	354af888c4	Create node-local DataCollection output folders Often times, compute nodes have local storage that is faster than the shared filesystem. Using node-local storage compared to the shared filesystem can also be advantageous to reduce the stress on such filesystem (which impacts all the users on a cluster). At the moment, `DataCollection::create_directory` creates the collection directory only on the global root rank (`myid == 0`) so that non-root nodes cannot write their per-rank ParaView and VisIt outputs when the path is not on the shared filesystem (e.g., on `/tmp` or `/scartch`). In this PR, I have the lowest rank on each shared-memory node (found via `MPI_COMM_TYPE_SHARED`) create the directory. When the filesystem is not shared, each node will have the folder where to write their outptu files. When the filesystem is shared, the extra mkdir() hits EEXIST, which is already tolerated, so behavior there is unchanged.	2026-07-22 16:05:59 -04: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