Helmholtz Multitrace novec is working.

.clang-format

@@ -0,0 +1,115 @@
+---
+Language:        Cpp
+# BasedOnStyle:  LLVM
+AccessModifierOffset: -2
+AlignAfterOpenBracket: Align
+AlignConsecutiveAssignments: false
+AlignConsecutiveDeclarations: false
+AlignEscapedNewlines: Right
+AlignOperands:   true
+AlignTrailingComments: true
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: false
+AllowShortCaseLabelsOnASingleLine: false
+AllowShortFunctionsOnASingleLine: All
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLoopsOnASingleLine: false
+AlwaysBreakAfterDefinitionReturnType: None
+AlwaysBreakAfterReturnType: None
+AlwaysBreakBeforeMultilineStrings: false
+AlwaysBreakTemplateDeclarations: false
+BinPackArguments: true
+BinPackParameters: true
+BraceWrapping:   
+  AfterClass:      false
+  AfterControlStatement: false
+  AfterEnum:       false
+  AfterFunction:   false
+  AfterNamespace:  false
+  AfterObjCDeclaration: false
+  AfterStruct:     false
+  AfterUnion:      false
+  AfterExternBlock: false
+  BeforeCatch:     false
+  BeforeElse:      false
+  IndentBraces:    false
+  SplitEmptyFunction: true
+  SplitEmptyRecord: true
+  SplitEmptyNamespace: true
+BreakBeforeBinaryOperators: None
+BreakBeforeBraces: Attach
+BreakBeforeInheritanceComma: false
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BreakConstructorInitializers: BeforeColon
+BreakAfterJavaFieldAnnotations: false
+BreakStringLiterals: true
+ColumnLimit:     80
+CommentPragmas:  '^ IWYU pragma:'
+CompactNamespaces: false
+ConstructorInitializerAllOnOneLineOrOnePerLine: false
+ConstructorInitializerIndentWidth: 4
+ContinuationIndentWidth: 4
+Cpp11BracedListStyle: true
+DerivePointerAlignment: false
+DisableFormat:   false
+ExperimentalAutoDetectBinPacking: false
+FixNamespaceComments: true
+ForEachMacros:   
+  - foreach
+  - Q_FOREACH
+  - BOOST_FOREACH
+IncludeBlocks:   Preserve
+IncludeCategories: 
+  - Regex:           '^"(llvm|llvm-c|clang|clang-c)/'
+    Priority:        2
+  - Regex:           '^(<|"(gtest|gmock|isl|json)/)'
+    Priority:        3
+  - Regex:           '.*'
+    Priority:        1
+IncludeIsMainRegex: '(Test)?$'
+IndentCaseLabels: false
+IndentPPDirectives: None
+IndentWidth:     2
+IndentWrappedFunctionNames: false
+JavaScriptQuotes: Leave
+JavaScriptWrapImports: true
+KeepEmptyLinesAtTheStartOfBlocks: true
+MacroBlockBegin: ''
+MacroBlockEnd:   ''
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: None
+ObjCBlockIndentWidth: 2
+ObjCSpaceAfterProperty: false
+ObjCSpaceBeforeProtocolList: true
+PenaltyBreakAssignment: 2
+PenaltyBreakBeforeFirstCallParameter: 19
+PenaltyBreakComment: 300
+PenaltyBreakFirstLessLess: 120
+PenaltyBreakString: 1000
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 60
+PointerAlignment: Right
+RawStringFormats: 
+  - Delimiter:       pb
+    Language:        TextProto
+    BasedOnStyle:    google
+ReflowComments:  true
+SortIncludes:    true
+SortUsingDeclarations: true
+SpaceAfterCStyleCast: false
+SpaceAfterTemplateKeyword: true
+SpaceBeforeAssignmentOperators: true
+SpaceBeforeParens: ControlStatements
+SpaceInEmptyParentheses: false
+SpacesBeforeTrailingComments: 1
+SpacesInAngles:  false
+SpacesInContainerLiterals: true
+SpacesInCStyleCastParentheses: false
+SpacesInParentheses: false
+SpacesInSquareBrackets: false
+Standard:        Cpp11
+TabWidth:        8
+UseTab:          Never
+...
+

bempp/api/operators/boundary/helmholtz.py

@@ -183,3 +183,89 @@ def hypersingular(
         device_interface,
         precision,
     )
+
+
+def multitrace_operator(
+    grid,
+    wavenumber,
+    parameters=None,
+    assembler="default_nonlocal",
+    device_interface=None,
+    precision=None,
+):
+    """Assemble the Helmholtz multitrace operator."""
+    from bempp.api.space import function_space
+    from bempp.api.operators import _add_wavenumber
+
+    if assembler != "multitrace_evaluator":
+        raise ValueError("Only multitrace evaluator supported.")
+
+    domain = function_space(grid, "P", 1)
+    range_ = domain
+    dual_to_range = domain
+
+    slp = single_layer(
+        domain,
+        range_,
+        dual_to_range,
+        wavenumber,
+        parameters,
+        "only_singular_part",
+        device_interface,
+        precision,
+    )
+
+    dlp = double_layer(
+        domain,
+        range_,
+        dual_to_range,
+        wavenumber,
+        parameters,
+        "only_singular_part",
+        device_interface,
+        precision,
+    )
+
+    adlp = adjoint_double_layer(
+        domain,
+        range_,
+        dual_to_range,
+        wavenumber,
+        parameters,
+        "only_singular_part",
+        device_interface,
+        precision,
+    )
+
+    hyp = hypersingular(
+        domain,
+        range_,
+        dual_to_range,
+        wavenumber,
+        parameters,
+        "only_singular_part",
+        device_interface,
+        precision,
+    )
+
+    options = {"COMPLEX_KERNEL": None}
+
+    _add_wavenumber(options, wavenumber)
+
+    singular_contribution = _np.array(
+        [[-dlp, slp], [hyp, adlp]], dtype=_np.object
+    )
+    return _common.create_multitrace_operator(
+        "helmholtz_multitrace",
+        [domain, domain],
+        [range_, range_],
+        [dual_to_range, dual_to_range],
+        parameters,
+        assembler,
+        options,
+        "helmholtz_multitrace",
+        singular_contribution,
+        device_interface,
+        precision,
+    )
+

bempp/api/operators/boundary/test/test_helmholtz.py

@@ -674,3 +674,128 @@ def test_helmholtz_hypersingular_evaluator(
 
     _np.testing.assert_allclose(actual, expected, rtol=tol)
 
+def test_helmholtz_hypersingular_complex_evaluator(
+    default_parameters, helpers, precision, device_interface
+):
+    """Test dense evaluator for the Helmholtz hypersingular op with complex wavenumber."""
+    from bempp.api import function_space
+    from bempp.api.operators.boundary.helmholtz import hypersingular
+
+    grid = helpers.load_grid("sphere")
+
+    space1 = function_space(grid, "P", 1)
+    space2 = function_space(grid, "P", 1)
+
+    discrete_op = hypersingular(
+        space1,
+        space1,
+        space2,
+        WAVENUMBER_COMPLEX,
+        assembler="dense_evaluator",
+        parameters=default_parameters,
+        device_interface=device_interface,
+        precision=precision,
+    ).assemble()
+
+    mat = hypersingular(
+        space1,
+        space1,
+        space2,
+        WAVENUMBER_COMPLEX,
+        assembler="dense",
+        parameters=default_parameters,
+        device_interface=device_interface,
+        precision=precision,
+    ).assemble()
+
+    x = _np.random.RandomState(0).randn(space1.global_dof_count)
+
+    actual = discrete_op @ x
+    expected = mat @ x
+
+    if precision == "single":
+        tol = 1e-4
+    else:
+        tol = 1e-12
+
+    _np.testing.assert_allclose(actual, expected, rtol=tol)
+
+
+def test_helmholtz_multitrace_sphere(default_parameters, helpers, device_interface, precision):
+    """Test Maxwell magnetic field on sphere."""
+    import bempp.api
+    from bempp.api import get_precision
+    from bempp.api import function_space
+    from bempp.api.shapes import regular_sphere
+    from bempp.api.operators.boundary.helmholtz import single_layer, double_layer, \
+            adjoint_double_layer, hypersingular
+    from bempp.api.assembly.blocked_operator import BlockedDiscreteOperator
+
+    # if precision == 'single':
+    #    pytest.skip("Test runs only in double precision mode.")
+
+    grid = helpers.load_grid('sphere')
+
+    op = bempp.api.operators.boundary.helmholtz.multitrace_operator(
+        grid,
+        2.5,
+        default_parameters,
+        assembler="multitrace_evaluator",
+        device_interface=device_interface,
+        precision=precision,
+    )
+
+    slp = single_layer(
+        op.domain_spaces[0],
+        op.range_spaces[0],
+        op.dual_to_range_spaces[0],
+        2.5,
+        parameters=default_parameters,
+        assembler="dense",
+        device_interface=device_interface,
+        precision=precision,
+    ).weak_form()
+
+    dlp = double_layer(
+        op.domain_spaces[0],
+        op.range_spaces[0],
+        op.dual_to_range_spaces[0],
+        2.5,
+        parameters=default_parameters,
+        assembler="dense",
+        device_interface=device_interface,
+        precision=precision,
+    ).weak_form()
+
+    adlp = adjoint_double_layer(
+        op.domain_spaces[0],
+        op.range_spaces[0],
+        op.dual_to_range_spaces[0],
+        2.5,
+        parameters=default_parameters,
+        assembler="dense",
+        device_interface=device_interface,
+        precision=precision,
+    ).weak_form()
+
+
+    hyp = hypersingular(
+        op.domain_spaces[0],
+        op.range_spaces[0],
+        op.dual_to_range_spaces[0],
+        2.5,
+        parameters=default_parameters,
+        assembler="dense",
+        device_interface=device_interface,
+        precision=precision,
+    ).weak_form()
+
+    expected = BlockedDiscreteOperator(_np.array([[-dlp, slp], [hyp, adlp]]))
+
+    rand = _np.random.RandomState(0)
+    x = rand.randn(expected.shape[1])
+
+    y_expected = expected @ x
+    y_actual = op.weak_form() @ x
+
+    _np.testing.assert_allclose(y_actual, y_expected, rtol=helpers.default_tolerance(precision))

bempp/core/dense_assembly_helpers.py

@@ -38,5 +38,7 @@ def choose_source_name_dense_multitrace_evaluator(compute_kernel):
 
     if compute_kernel == "maxwell_multitrace":
         return "evaluate_dense_vector_maxwell_multitrace"
+    if compute_kernel == "helmholtz_multitrace":
+        return "evaluate_dense_vector_helmholtz_multitrace"
 
     raise ValueError("Unknown compute kernel identifier.")

bempp/core/sources/include/bempp_base_types.h

@@ -57,7 +57,10 @@ typedef struct Geometry
     REALTYPE volume;
 } Geometry;
 
+/* Access an element of a vector type. */
 #define VEC_ELEMENT(A, INDEX) ((REALTYPE*)&A)[INDEX]
 
+/* Print out a complex variable for debugging. */
+#define PRINT_COMPLEX(A, INFO) printf(INFO" %e %e\n", A[0], A[1])
 
 #endif

bempp/core/sources/kernels/evaluate_dense_vector_helmholtz_hypersingular_regular_novec.cl

@@ -51,6 +51,7 @@ __kernel void evaluate_dense_helmholtz_hypersingular_novec(
 
   REALTYPE testInv[2][2];
   REALTYPE trialInv[2][2];
+  REALTYPE tmp[2];
 
   REALTYPE3 trialCurl[3];
   REALTYPE3 testCurl[3];
@@ -245,14 +246,16 @@ __local REALTYPE localResult[WORKGROUP_SIZE][3][3];
 
   for (i = 0; i < 3; ++i)
     for (j = 0; j < 3; ++j) {
+      tmp[0] = shapeIntegral[i][j][0];
+      tmp[1] = shapeIntegral[i][j][1];
       shapeIntegral[i][j][0] =
-          -(wavenumberProduct[0] * shapeIntegral[i][j][0] -
-            wavenumberProduct[1] * shapeIntegral[i][j][1]) *
+          -(wavenumberProduct[0] * tmp[0] -
+            wavenumberProduct[1] * tmp[1]) *
               normalProduct +
           firstTermIntegral[0] * basisProduct[i][j];
       shapeIntegral[i][j][1] =
-          -(wavenumberProduct[0] * shapeIntegral[i][j][1] +
-            wavenumberProduct[1] * shapeIntegral[i][j][0]) *
+          -(wavenumberProduct[0] * tmp[1] +
+            wavenumberProduct[1] * tmp[0]) *
               normalProduct +
           firstTermIntegral[1] * basisProduct[i][j];
     }