C++ BoundaryOperator类代码示例

BOOST_AUTO_TEST_CASE_TEMPLATE(symmetric_matches_nonsymmetric_in_aca_mode,
                              ValueType, result_types)
{
    typedef ValueType RT;
    typedef typename Fiber::ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    if (boost::is_same<RT, std::complex<float> >::value) {
        // The AHMED support for single-precision complex symmetric matrices
        // is broken
        BOOST_CHECK(true);
        return;
    }

    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(
                                params, "../../examples/meshes/sphere-h-0.4.msh",
                                false /* verbose */);

    PiecewiseLinearContinuousScalarSpace<BFT> pwiseLinears(grid);
    PiecewiseConstantScalarSpace<BFT> pwiseConstants(grid);

    AssemblyOptions assemblyOptions;
    assemblyOptions.setVerbosityLevel(VerbosityLevel::LOW);
    AcaOptions acaOptions;
    acaOptions.minimumBlockSize = 4;
    assemblyOptions.switchToAcaMode(acaOptions);
    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setRelativeQuadratureOrder(4);
    accuracyOptions.doubleSingular.setRelativeQuadratureOrder(2);
    NumericalQuadratureStrategy<BFT, RT> quadStrategy(accuracyOptions);

    Context<BFT, RT> context(make_shared_from_ref(quadStrategy), assemblyOptions);

    const RT waveNumber = initWaveNumber<RT>();

    BoundaryOperator<BFT, RT> opNonsymmetric =
        modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, RT, RT>(
            make_shared_from_ref(context),
            make_shared_from_ref(pwiseLinears),
            make_shared_from_ref(pwiseConstants),
            make_shared_from_ref(pwiseLinears),
            waveNumber,
            "", NO_SYMMETRY);
    BoundaryOperator<BFT, RT> opSymmetric =
        modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, RT, RT>(
            make_shared_from_ref(context),
            make_shared_from_ref(pwiseLinears),
            make_shared_from_ref(pwiseConstants),
            make_shared_from_ref(pwiseLinears),
            waveNumber,
            "", SYMMETRIC);

    arma::Mat<RT> matNonsymmetric = opNonsymmetric.weakForm()->asMatrix();
    arma::Mat<RT> matSymmetric = opSymmetric.weakForm()->asMatrix();

    BOOST_CHECK(check_arrays_are_close<RT>(
                    matNonsymmetric, matSymmetric, 2 * acaOptions.eps));
}

BOOST_AUTO_TEST_CASE_TEMPLATE(aca_of_synthetic_maxwell_single_layer_operator_agrees_with_dense_assembly_in_asymmetric_case,
                              ValueType, complex_result_types)
{
    typedef ValueType RT;
    typedef typename ScalarTraits<ValueType>::RealType RealType;
    typedef RealType BFT;

    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(
        params, "../../meshes/sphere-h-0.4.msh", false /* verbose */);

    RT waveNumber = initWaveNumber<RT>();

    shared_ptr<Space<BFT> > vectorPwiseLinears(
        new RaviartThomas0VectorSpace<BFT>(grid));
    shared_ptr<Space<BFT> > vectorPwiseLinears2(
        new RaviartThomas0VectorSpace<BFT>(grid));

    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setRelativeQuadratureOrder(2);
    accuracyOptions.singleRegular.setRelativeQuadratureOrder(2);
    shared_ptr<NumericalQuadratureStrategy<BFT, RT> > quadStrategy(
                new NumericalQuadratureStrategy<BFT, RT>(accuracyOptions));

    AssemblyOptions assemblyOptionsDense;
    assemblyOptionsDense.setVerbosityLevel(VerbosityLevel::LOW);
    shared_ptr<Context<BFT, RT> > contextDense(
        new Context<BFT, RT>(quadStrategy, assemblyOptionsDense));

    BoundaryOperator<BFT, RT> opDense =
            maxwell3dSingleLayerBoundaryOperator<BFT>(
                contextDense,
                vectorPwiseLinears, vectorPwiseLinears, vectorPwiseLinears,
                waveNumber);
    arma::Mat<RT> weakFormDense = opDense.weakForm()->asMatrix();

    AssemblyOptions assemblyOptionsAca;
    assemblyOptionsAca.setVerbosityLevel(VerbosityLevel::LOW);
    AcaOptions acaOptions;
    acaOptions.mode = AcaOptions::LOCAL_ASSEMBLY;
    assemblyOptionsAca.switchToAcaMode(acaOptions);
    shared_ptr<Context<BFT, RT> > contextAca(
        new Context<BFT, RT>(quadStrategy, assemblyOptionsAca));

    // Internal domain different from dualToRange
    BoundaryOperator<BFT, RT> opAca =
            maxwell3dSingleLayerBoundaryOperator<BFT>(
                contextAca,
                vectorPwiseLinears, vectorPwiseLinears, vectorPwiseLinears2,
                waveNumber);
    arma::Mat<RT> weakFormAca = opAca.weakForm()->asMatrix();

    BOOST_CHECK(check_arrays_are_close<ValueType>(
                    weakFormDense, weakFormAca, 2. * acaOptions.eps));
}

BOOST_AUTO_TEST_CASE_TEMPLATE(interpolated_matches_noniterpolated,
                              BasisFunctionType, basis_function_types)
{
    typedef BasisFunctionType BFT;
    typedef typename Fiber::ScalarTraits<BFT>::ComplexType RT;
    typedef typename Fiber::ScalarTraits<BFT>::ComplexType KT;
    typedef typename Fiber::ScalarTraits<BFT>::RealType CT;

    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(
                params, "meshes/two_disjoint_triangles.msh",
                false /* verbose */);

    PiecewiseLinearContinuousScalarSpace<BFT> pwiseLinears(grid);
    PiecewiseConstantScalarSpace<BFT> pwiseConstants(grid);

    AssemblyOptions assemblyOptions;
    assemblyOptions.setVerbosityLevel(VerbosityLevel::LOW);
    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setAbsoluteQuadratureOrder(5);
    accuracyOptions.doubleSingular.setAbsoluteQuadratureOrder(5);
    NumericalQuadratureStrategy<BFT, RT> quadStrategy(accuracyOptions);

    Context<BFT, RT> context(make_shared_from_ref(quadStrategy), assemblyOptions);

    const KT waveNumber(3.23, 0.31);

    BoundaryOperator<BFT, RT> opNoninterpolated =
            modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, KT, RT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                waveNumber,
                "", NO_SYMMETRY,
                false);
    BoundaryOperator<BFT, RT> opInterpolated =
            modifiedHelmholtz3dSingleLayerBoundaryOperator<BFT, KT, RT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                waveNumber,
                "", NO_SYMMETRY,
                true);

    arma::Mat<RT> matNoninterpolated = opNoninterpolated.weakForm()->asMatrix();
    arma::Mat<RT> matInterpolated = opInterpolated.weakForm()->asMatrix();

    const CT eps = std::numeric_limits<CT>::epsilon();
    BOOST_CHECK(check_arrays_are_close<RT>(
                    matNoninterpolated, matInterpolated, 100 * eps));
}

void BlockedOperatorStructure<BasisFunctionType, ResultType>::setBlock(
    size_t row, size_t column,
    const BoundaryOperator<BasisFunctionType, ResultType> &op) {
  if (!op.isInitialized())
    resetBlock(row, column);
  else
    m_blocks[Key(row, column)] = op;
}

BoundaryOperator<BasisFunctionType, ResultType>
pseudoinverse(const BoundaryOperator<BasisFunctionType, ResultType>& boundaryOp)
{
    typedef AbstractBoundaryOperatorPseudoinverse<BasisFunctionType, ResultType>
            Pinv;
    return BoundaryOperator<BasisFunctionType, ResultType>(
                boundaryOp.context(),
                boost::make_shared<Pinv>(boundaryOp));
}

AbstractBoundaryOperatorPseudoinverse<BasisFunctionType, ResultType>::
AbstractBoundaryOperatorPseudoinverse(
        // TODO: add a solver argument specifying how to calculate the pseudoinv.
        const BoundaryOperator<BasisFunctionType, ResultType>& operatorToInvert,
        const shared_ptr<const Space<BasisFunctionType> >& dualToRange) :
    Base(operatorToInvert.range(), operatorToInvert.domain(),dualToRange,
         "pinv(" + operatorToInvert.label() + ")",
         throwIfUninitialized(operatorToInvert,
                              "AbstractBoundaryOperatorPseudoinverse::"
                              "AbstractBoundaryOperatorPseudoinverse(): "
                              "the boundary operator to be inverted must be "
                              "initialized"
                              ).abstractOperator()->symmetry()),
    m_operator(operatorToInvert),
    m_id(boost::make_shared<AbstractBoundaryOperatorPseudoinverseId<
         BasisFunctionType, ResultType> >(
             operatorToInvert))
{
}

AdjointAbstractBoundaryOperator<BasisFunctionType, ResultType>::
AdjointAbstractBoundaryOperator(
        const BoundaryOperator<BasisFunctionType, ResultType>& boundaryOp,
        int symmetry) :
    Base(boundaryOp.dualToRange(),
         boundaryOp.domain() == boundaryOp.range() ?
             boundaryOp.dualToRange() :
             // assume that domain == dualToRange, we'll verify it
             // in the body of the constructor
             boundaryOp.range(),
         boundaryOp.domain(),
         "adj(" + boundaryOp.label() + ")",
         symmetry & AUTO_SYMMETRY ?
             boundaryOp.abstractOperator()->symmetry() :
             symmetry),
    m_operator(boundaryOp)
{
    if (boost::is_complex<BasisFunctionType>())
        throw std::logic_error(
            "AdjointAbstractBoundaryOperator(): Taking the adjoint of operators "
            "acting on complex-valued basis functions is not supported");
    if (boundaryOp.domain() != boundaryOp.range() &&
            boundaryOp.domain() != boundaryOp.dualToRange())
        throw std::runtime_error(
                "AdjointAbstractBoundaryOperator::"
                "AdjointAbstractBoundaryOperator(): "
                "Dual to the domain of the operator to invert cannot be determined "
                "since the domain is different from both "
                "the range and the space dual to its range");
}

AbstractBoundaryOperatorPseudoinverseId<BasisFunctionType, ResultType>::
AbstractBoundaryOperatorPseudoinverseId(
        const BoundaryOperator<BasisFunctionType, ResultType>& operatorToInvert) :
    m_operatorToInvertId(operatorToInvert.abstractOperator()->id())
{
}

BOOST_AUTO_TEST_CASE_TEMPLATE(works, BasisFunctionType, basis_function_types)
{
    typedef BasisFunctionType BFT;
    typedef typename Fiber::ScalarTraits<BFT>::ComplexType RT;
    typedef typename Fiber::ScalarTraits<BFT>::RealType CT;
    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(
                params, "meshes/two_disjoint_triangles.msh",
                false /* verbose */);

    PiecewiseLinearContinuousScalarSpace<BFT> pwiseLinears(grid);
    PiecewiseConstantScalarSpace<BFT> pwiseConstants(grid);

    AssemblyOptions assemblyOptions;
    assemblyOptions.setVerbosityLevel(VerbosityLevel::LOW);
    AccuracyOptions accuracyOptions;
    accuracyOptions.doubleRegular.setAbsoluteQuadratureOrder(5);
    accuracyOptions.doubleSingular.setAbsoluteQuadratureOrder(5);
    NumericalQuadratureStrategy<BFT, RT> quadStrategy(accuracyOptions);

    Context<BFT, RT> context(make_shared_from_ref(quadStrategy), assemblyOptions);

    const RT waveNumber(1.23, 0.31);

    BoundaryOperator<BFT, RT> slpOpConstants = Bempp::helmholtz3dSingleLayerBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseConstants),
                make_shared_from_ref(pwiseConstants),
                make_shared_from_ref(pwiseConstants),
                waveNumber);
    BoundaryOperator<BFT, RT> slpOpLinears = helmholtz3dSingleLayerBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                waveNumber);
    BoundaryOperator<BFT, RT> hypOp = helmholtz3dHypersingularBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                make_shared_from_ref(pwiseLinears),
                waveNumber);

    // Get the matrix repr. of the hypersingular operator
    arma::Mat<RT> hypMat = hypOp.weakForm()->asMatrix();

    // Construct the expected hypersingular operator matrix. For this, we need:

    // * the surface curls of all basis functions (which are constant)
    typedef Fiber::SurfaceCurl3dElementaryFunctor<CT> ElementaryFunctor;
    typedef Fiber::ElementaryBasisTransformationFunctorWrapper<ElementaryFunctor> Functor;
    Functor functor;
    size_t basisDeps = 0, geomDeps = 0;
    functor.addDependencies(basisDeps, geomDeps);

    arma::Mat<CT> points(2, 1);
    points.fill(0.);

    typedef Fiber::PiecewiseLinearContinuousScalarBasis<3, BFT> Basis;
    Basis basis;
    Fiber::BasisData<BFT> basisData;
    basis.evaluate(basisDeps, points, Fiber::ALL_DOFS, basisData);

    Fiber::GeometricalData<CT> geomData[2];
    std::unique_ptr<GridView> view = grid->leafView();
    std::unique_ptr<EntityIterator<0> > it = view->entityIterator<0>();
    it->entity().geometry().getData(geomDeps, points, geomData[0]);
    it->next();
    it->entity().geometry().getData(geomDeps, points, geomData[1]);

    Fiber::DefaultCollectionOfBasisTransformations<Functor> transformations(functor);

    Fiber::CollectionOf3dArrays<BFT> surfaceCurl[2];
    transformations.evaluate(basisData, geomData[0], surfaceCurl[0]);
    transformations.evaluate(basisData, geomData[1], surfaceCurl[1]);

    // * the single-layer potential matrix for constant basis functions
    arma::Mat<RT> slpMatConstants = slpOpConstants.weakForm()->asMatrix();
    // * the single-layer potential matrix for linear basis functions
    arma::Mat<RT> slpMatLinears = slpOpLinears.weakForm()->asMatrix();

    arma::Mat<RT> expectedHypMat(6, 6);
    for (size_t testElement = 0; testElement < 2; ++testElement)
        for (size_t trialElement = 0; trialElement < 2; ++trialElement)
            for (size_t r = 0; r < 3; ++r)
                for (size_t c = 0; c < 3; ++c) {
                    RT curlMultiplier = 0.;
                    for (size_t dim = 0; dim < 3; ++dim)
                        curlMultiplier += surfaceCurl[testElement][0](dim, r, 0) *
                                surfaceCurl[trialElement][0](dim, c, 0);
                    RT valueMultiplier = 0.;
                    for (size_t dim = 0; dim < 3; ++dim)
                        valueMultiplier += geomData[testElement].normals(dim, 0) *
                            geomData[trialElement].normals(dim, 0);
                    expectedHypMat(3 * testElement + r, 3 * trialElement + c) =
                            curlMultiplier * slpMatConstants(testElement, trialElement) -
                            waveNumber * waveNumber * valueMultiplier *
                            slpMatLinears(3 * testElement + r, 3 * trialElement + c);
                }
//.........这里部分代码省略.........

int main(int argc, char* argv[])
{
    // Process command-line args

    if (argc < 7 || argc % 2 != 1) {
        std::cout << "Solve a Maxwell Dirichlet problem in an exterior domain.\n"
                  << "Usage: " << argv[0]
                  << " <mesh_file> <n_threads> <aca_eps> <solver_tol>"
                  << " <singular_order_increment>"
                  << " [<regular_order_increment_1> <min_relative_distance_1>]"
                  << " [<regular_order_increment_2> <min_relative_distance_2>]"
                  << " [...] <regular_order_increment_n>"
                  << std::endl;
        return 1;
    }
    int maxThreadCount = atoi(argv[2]);
    double acaEps = atof(argv[3]);
    double solverTol = atof(argv[4]);
    int singOrderIncrement = atoi(argv[5]);
    if (acaEps > 1. || acaEps < 0.) {
        std::cout << "Invalid aca_eps: " << acaEps << std::endl;
        return 1;
    }
    if (solverTol > 1. || solverTol < 0.) {
        std::cout << "Invalid solver_tol: " << solverTol << std::endl;
        return 1;
    }

    AccuracyOptionsEx accuracyOptions;
    std::vector<double> maxNormDists;
    std::vector<int> orderIncrements;
    for (int i = 6; i < argc - 1; i += 2) {
        orderIncrements.push_back(atoi(argv[i]));
        maxNormDists.push_back(atof(argv[i + 1]));
    }
    orderIncrements.push_back(atoi(argv[argc - 1]));
    accuracyOptions.setDoubleRegular(maxNormDists, orderIncrements);
    accuracyOptions.setDoubleSingular(singOrderIncrement);
    accuracyOptions.setSingleRegular(2);

    // Load mesh

    GridParameters params;
    params.topology = GridParameters::TRIANGULAR;
    shared_ptr<Grid> grid = GridFactory::importGmshGrid(params, argv[1]);

    // Initialize the space

    RaviartThomas0VectorSpace<BFT> HdivSpace(grid);

    // Set assembly mode and options

    AssemblyOptions assemblyOptions;
    assemblyOptions.setMaxThreadCount(maxThreadCount);
    if (acaEps > 0.) {
        AcaOptions acaOptions;
        acaOptions.eps = acaEps;
        assemblyOptions.switchToAcaMode(acaOptions);
    }

    NumericalQuadratureStrategy<BFT, RT> quadStrategy(accuracyOptions);
    Context<BFT, RT> context(make_shared_from_ref(quadStrategy), assemblyOptions);

    // Construct operators

    BoundaryOperator<BFT, RT> slpOp = maxwell3dSingleLayerBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                k,
                "SLP");
    BoundaryOperator<BFT, RT> dlpOp = maxwell3dDoubleLayerBoundaryOperator<BFT>(
                make_shared_from_ref(context),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                k,
                "DLP");
    BoundaryOperator<BFT, RT> idOp = maxwell3dIdentityOperator<BFT, RT>(
                make_shared_from_ref(context),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                "Id");

    // Construct a grid function representing the Dirichlet data

    GridFunction<BFT, RT> dirichletData(
                make_shared_from_ref(context),
                make_shared_from_ref(HdivSpace),
                make_shared_from_ref(HdivSpace),
                surfaceNormalDependentFunction(DirichletData()));

    dirichletData.exportToVtk(VtkWriter::CELL_DATA, "Dirichlet_data",
                              "input_dirichlet_data_cell");
    dirichletData.exportToVtk(VtkWriter::VERTEX_DATA, "Dirichlet_data",
                              "input_dirichlet_data_vertex");

    // Construct a grid function representing the right-hand side
//.........这里部分代码省略.........

BoundaryOperator<BasisFunctionType, ResultType>
modifiedHelmholtz3dSyntheticHypersingularBoundaryOperator(
    const shared_ptr<const Context<BasisFunctionType, ResultType>> &context,
    const shared_ptr<const Space<BasisFunctionType>> &domain,
    const shared_ptr<const Space<BasisFunctionType>> &range,
    const shared_ptr<const Space<BasisFunctionType>> &dualToRange,
    KernelType waveNumber, std::string label, int internalSymmetry,
    bool useInterpolation, int interpPtsPerWavelength,
    const BoundaryOperator<BasisFunctionType, ResultType> &externalSlp) {
  typedef typename ScalarTraits<BasisFunctionType>::RealType CoordinateType;

  typedef Fiber::ScalarFunctionValueFunctor<CoordinateType> ValueFunctor;
  typedef Fiber::ScalarFunctionValueTimesNormalFunctor<CoordinateType>
      ValueTimesNormalFunctor;
  typedef Fiber::SurfaceCurl3dFunctor<CoordinateType> CurlFunctor;
  typedef Fiber::SingleComponentTestTrialIntegrandFunctor<
      BasisFunctionType, ResultType> IntegrandFunctor;

  typedef GeneralElementaryLocalOperator<BasisFunctionType, ResultType> LocalOp;
  typedef SyntheticIntegralOperator<BasisFunctionType, ResultType> SyntheticOp;

  if (!domain || !range || !dualToRange)
    throw std::invalid_argument(
        "modifiedHelmholtz3dSyntheticHypersingularBoundaryOperator(): "
        "domain, range and dualToRange must not be null");

  shared_ptr<const Space<BasisFunctionType>> newDomain = domain;
  shared_ptr<const Space<BasisFunctionType>> newDualToRange = dualToRange;

  bool isBarycentric =
      (domain->isBarycentric() || dualToRange->isBarycentric());

  if (isBarycentric) {
    newDomain = domain->barycentricSpace(domain);
    newDualToRange = dualToRange->barycentricSpace(dualToRange);
  }

  shared_ptr<const Context<BasisFunctionType, ResultType>> internalContext,
      auxContext;
  SyntheticOp::getContextsForInternalAndAuxiliaryOperators(
      context, internalContext, auxContext);
  shared_ptr<const Space<BasisFunctionType>> internalTrialSpace =
      newDomain->discontinuousSpace(newDomain);
  shared_ptr<const Space<BasisFunctionType>> internalTestSpace =
      newDualToRange->discontinuousSpace(newDualToRange);

  // Note: we don't really need to care about ranges and duals to domains of
  // the internal operator. The only range space that matters is that of the
  // leftmost operator in the product.

  const char xyz[] = "xyz";
  const size_t dimWorld = 3;

  if (label.empty())
    label =
        AbstractBoundaryOperator<BasisFunctionType, ResultType>::uniqueLabel();

  BoundaryOperator<BasisFunctionType, ResultType> slp;
  if (!externalSlp.isInitialized()) {

    slp = modifiedHelmholtz3dSingleLayerBoundaryOperator<
        BasisFunctionType, KernelType, ResultType>(
        internalContext, internalTrialSpace,
        internalTestSpace /* or whatever */, internalTestSpace, waveNumber,
        "(" + label + ")_internal_SLP", internalSymmetry, useInterpolation,
        interpPtsPerWavelength);
  } else {

    slp = externalSlp;
  }

  // symmetry of the decomposition
  int syntheseSymmetry = 0; // symmetry of the decomposition
  if (newDomain == newDualToRange && internalTrialSpace == internalTestSpace)
    syntheseSymmetry =
        HERMITIAN | (boost::is_complex<BasisFunctionType>() ? 0 : SYMMETRIC);

  std::vector<BoundaryOperator<BasisFunctionType, ResultType>> testLocalOps;
  std::vector<BoundaryOperator<BasisFunctionType, ResultType>> trialLocalOps;
  testLocalOps.resize(dimWorld);
  for (size_t i = 0; i < dimWorld; ++i)
    testLocalOps[i] = BoundaryOperator<BasisFunctionType, ResultType>(
        auxContext, boost::make_shared<LocalOp>(
                        internalTestSpace, range, newDualToRange,
                        ("(" + label + ")_test_curl_") + xyz[i], NO_SYMMETRY,
                        CurlFunctor(), ValueFunctor(), IntegrandFunctor(i, 0)));

  if (!syntheseSymmetry) {
    trialLocalOps.resize(dimWorld);
    for (size_t i = 0; i < dimWorld; ++i)
      trialLocalOps[i] = BoundaryOperator<BasisFunctionType, ResultType>(
          auxContext,
          boost::make_shared<LocalOp>(
              newDomain, internalTrialSpace /* or whatever */,
              internalTrialSpace, ("(" + label + ")_trial_curl_") + xyz[i],
              NO_SYMMETRY, ValueFunctor(), CurlFunctor(),
              IntegrandFunctor(0, i)));
  }
  // It might be more prudent to distinguish between the symmetry of the total
  // operator and the symmetry of the decomposition
//.........这里部分代码省略.........

BoundaryOperator<BasisFunctionType, ResultType>
modifiedHelmholtz3dHypersingularBoundaryOperator(
    const shared_ptr<const Context<BasisFunctionType, ResultType>> &context,
    const shared_ptr<const Space<BasisFunctionType>> &domain,
    const shared_ptr<const Space<BasisFunctionType>> &range,
    const shared_ptr<const Space<BasisFunctionType>> &dualToRange,
    KernelType waveNumber, const std::string &label, int symmetry,
    bool useInterpolation, int interpPtsPerWavelength,
    const BoundaryOperator<BasisFunctionType, ResultType> &externalSlp) {
  const AssemblyOptions &assemblyOptions = context->assemblyOptions();
  if ((assemblyOptions.assemblyMode() == AssemblyOptions::ACA &&
       assemblyOptions.acaOptions().mode == AcaOptions::LOCAL_ASSEMBLY) ||
      externalSlp.isInitialized())
    return modifiedHelmholtz3dSyntheticHypersingularBoundaryOperator(
        context, domain, range, dualToRange, waveNumber, label, symmetry,
        useInterpolation, interpPtsPerWavelength, externalSlp);

  typedef typename ScalarTraits<BasisFunctionType>::RealType CoordinateType;

  typedef Fiber::ModifiedHelmholtz3dHypersingularKernelFunctor<KernelType>
      NoninterpolatedKernelFunctor;
  typedef Fiber::ModifiedHelmholtz3dHypersingularKernelInterpolatedFunctor<
      KernelType> InterpolatedKernelFunctor;
  typedef Fiber::ModifiedHelmholtz3dHypersingularTransformationFunctor<
      CoordinateType> TransformationFunctor;
  typedef Fiber::ModifiedHelmholtz3dHypersingularIntegrandFunctor2<
      BasisFunctionType, KernelType, ResultType> IntegrandFunctor;

  typedef Fiber::ModifiedHelmholtz3dHypersingularTransformationFunctor2<
      CoordinateType> TransformationFunctorWithBlas;

  typedef Fiber::
      ModifiedHelmholtz3dHypersingularOffDiagonalInterpolatedKernelFunctor<
          KernelType> OffDiagonalInterpolatedKernelFunctor;
  typedef Fiber::ModifiedHelmholtz3dHypersingularOffDiagonalKernelFunctor<
      KernelType> OffDiagonalNoninterpolatedKernelFunctor;
  typedef Fiber::ScalarFunctionValueFunctor<CoordinateType>
      OffDiagonalTransformationFunctor;
  typedef Fiber::SimpleTestScalarKernelTrialIntegrandFunctorExt<
      BasisFunctionType, KernelType, ResultType, 1> OffDiagonalIntegrandFunctor;

  CoordinateType maxDistance_ =
      static_cast<CoordinateType>(1.1) *
      maxDistance(*domain->grid(), *dualToRange->grid());

  shared_ptr<Fiber::TestKernelTrialIntegral<BasisFunctionType, KernelType,
                                            ResultType>> integral,
      offDiagonalIntegral;
  if (shouldUseBlasInQuadrature(assemblyOptions, *domain, *dualToRange)) {
    integral.reset(new Fiber::TypicalTestScalarKernelTrialIntegral<
        BasisFunctionType, KernelType, ResultType>());
    offDiagonalIntegral = integral;
  } else {
    integral.reset(new Fiber::DefaultTestKernelTrialIntegral<IntegrandFunctor>(
        IntegrandFunctor()));
    offDiagonalIntegral.reset(
        new Fiber::DefaultTestKernelTrialIntegral<OffDiagonalIntegrandFunctor>(
            OffDiagonalIntegrandFunctor()));
  }

  typedef GeneralHypersingularIntegralOperator<BasisFunctionType, KernelType,
                                               ResultType> Op;
  shared_ptr<Op> newOp;
  if (shouldUseBlasInQuadrature(assemblyOptions, *domain, *dualToRange)) {
    shared_ptr<Fiber::TestKernelTrialIntegral<BasisFunctionType, KernelType,
                                              ResultType>> integral,
        offDiagonalIntegral;
    integral.reset(new Fiber::TypicalTestScalarKernelTrialIntegral<
        BasisFunctionType, KernelType, ResultType>());
    offDiagonalIntegral = integral;
    if (useInterpolation)
      newOp.reset(new Op(
          domain, range, dualToRange, label, symmetry,
          InterpolatedKernelFunctor(waveNumber, maxDistance_,
                                    interpPtsPerWavelength),
          TransformationFunctorWithBlas(), TransformationFunctorWithBlas(),
          integral, OffDiagonalInterpolatedKernelFunctor(
                        waveNumber, maxDistance_, interpPtsPerWavelength),
          OffDiagonalTransformationFunctor(),
          OffDiagonalTransformationFunctor(), offDiagonalIntegral));
    else
      newOp.reset(new Op(
          domain, range, dualToRange, label, symmetry,
          NoninterpolatedKernelFunctor(waveNumber),
          TransformationFunctorWithBlas(), TransformationFunctorWithBlas(),
          integral, OffDiagonalNoninterpolatedKernelFunctor(waveNumber),
          OffDiagonalTransformationFunctor(),
          OffDiagonalTransformationFunctor(), offDiagonalIntegral));
  } else { // no blas
    if (useInterpolation)
      newOp.reset(new Op(
          domain, range, dualToRange, label, symmetry,
          InterpolatedKernelFunctor(waveNumber, maxDistance_,
                                    interpPtsPerWavelength),
          TransformationFunctor(), TransformationFunctor(), IntegrandFunctor(),
          OffDiagonalInterpolatedKernelFunctor(waveNumber, maxDistance_,
                                               interpPtsPerWavelength),
          OffDiagonalTransformationFunctor(),
          OffDiagonalTransformationFunctor(), OffDiagonalIntegrandFunctor()));
    else
//.........这里部分代码省略.........