C++ BoundaryOperator::weakForm方法代码示例

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

BOOST_AUTO_TEST_CASE_TEMPLATE(aca_of_synthetic_modified_helmholtz_hypersingular_operator_agrees_with_dense_assembly_in_symmetric_case,
                              ValueType, 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> > pwiseConstants(
        new PiecewiseConstantScalarSpace<BFT>(grid));
    shared_ptr<Space<BFT> > pwiseLinears(
        new PiecewiseLinearContinuousScalarSpace<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 =
            modifiedHelmholtz3dHypersingularBoundaryOperator<BFT, RT, RT>(
                contextDense, pwiseLinears, pwiseConstants, pwiseLinears,
                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));

    BoundaryOperator<BFT, RT> opAca =
            modifiedHelmholtz3dHypersingularBoundaryOperator<BFT, RT, RT>(
                contextAca, pwiseLinears, pwiseConstants, pwiseLinears,
                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(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
//.........这里部分代码省略.........